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Abstract 

The  Kolmogorov,  Arnold,  and  Moser  (KAM)  Theory  was  developed  in  the  1960s 
but  only  in  the  last  decade  has  it  been  applied  to  Earth  orbiting  satellites.  Physical 
state  variables  of  position  and  velocity  are  transformed  into  KAM  Torus  variables. 
The  KAM  Torus  is  a  geometrical  structure  similar  to  that  of  a  multi-dimensional 
donut.  The  satellite’s  motion  can  be  described  as  traversing  the  surface  of  this  donut. 
There  are  two  primary  advantages  of  this  transformation:  1)  The  new  generalized 
coordinates  which  are  analogous  with  mean  anomaly,  right  ascension  of  the  ascending 
node,  and  argument  of  perigee,  increment  linearly  with  time,  and  2)  Perturbations 
due  to  the  Earth’s  geopotential  are  already  embedded  in  a  given  torus  to  an  arbitrary 
geopotential  order.  This  study  examines  methods  to  describe  perturbed  satellite  mo¬ 
tion  near  a  reference  KAM  Torus.  The  perturbations  addressed  in  this  thesis  are 
atmospheric  air  drag  and  third-body  effects  from  the  Moon. 

Perturbed  motion  was  integrated  and  compared  against  the  unperturbed  refer¬ 
ence  torus  motion.  For  a  sample  orbit,  expressions  were  numerically  derived  that  allow 
the  modification  of  the  reference  torus  to  allow  prediction  of  the  perturbed  motion 
due  to  drag.  It  was  shown  that  in  the  case  of  third-body  lunar  effects,  the  differences 
between  perturbed  motion  and  the  reference  torus  motion  cannot  be  generalized.  In¬ 
stead,  there  seems  to  be  both  evidence  and  motivation  behind  attempting  to  embed 
the  lunar  dynamics  in  the  construction  of  Earth  satellite  KAM  Tori. 
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Effects  of  Air  Drag  and  Lunar  Third-Body 
Pertubrations  on  Orbital  Motion  Near 
a  Reference  KAM  Torus 

I.  Introduction 

Though  ancient  peoples  have  been  watching  the  nighttime  sky  for  millennia,  the 
study  of  orbital  mechanics  as  we  know  it  began  with  such  giants  as  Nicholas  Coper¬ 
nicus,  Johannes  Kepler,  and  Galileo  Galilei  in  the  17th  century.  These  great  thinkers 
were  simply  trying  to  reconcile  the  motion  of  the  heavens  with  rudimentary  thoughts 
such  as  the  earth  being  flat  and  the  center  of  the  known  universe.  Their  observations 
reasonably  focused  on  heavenly  bodies  such  as  the  Sun,  Moon,  and  planets. 

Flip  a  few  pages  forward  in  the  calendar  to  the  twentieth  century,  and  humanity 
can  be  found  leveraging  orbital  mechanics  for  amazing  applications.  The  Soviet  Union 
placed  the  first  artificial  satellite  in  orbit  around  the  Earth  in  1957.  The  United 
States  of  America  landed  men  on  the  moon  in  1969.  Again,  fast  forward  in  time,  but 
now  to  the  current  day.  The  nations  of  the  world  utilize  space  both  commercially 
and  militarily.  Commercial  applications  include  Sun  and  Earth  weather  observation, 
communication,  entertainment,  and  imaging.  Military  applications  include  various 
forms  of  intelligence  gathering,  communications,  precision  navigation  and  timing,  and 
space  situational  awareness. 

1 . 1  Motivation 

Somewhere  in  the  transition  from  predicting  the  motion  of  heavenly  bodies  to 
the  placement  of  artificial  satellites  in  low-earth  orbit,  it  became  necessary  to  consider 
how  classical  satellite  orbits  as  defined  by  Kepler  are  perturbed  by  small  forces  over 
time.  Current  day  space  applications  require  high  accuracy  knowledge  of  a  spacecraft’s 
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position  as  well  as  the  ability  to  predict  the  spacecraft’s  state  at  some  time  in  the 
future.  Today,  this  is  achieved  by  collecting  observations  of  position  and  velocity  of 
the  satellite  using  ground-based  sensors.  These  sensors  include  optical  sensors  and 
the  “Space  Fence,”  a  very  high  frequency  radar  network,  all  owned  and  operated  by 
the  United  States.  The  observations  are  combined  to  produce  a  description,  or  state, 
of  the  spacecraft’s  orbit,  called  orbital  elements.  The  spacecraft’s  state  can  then  be 
propagated  to  some  arbitrary  time  in  the  future  by  means  of  numerical  integration. 

Unfortunately,  as  the  “space  is  big”  theory  for  on-orbit  disposal  of  dead  satel¬ 
lites,  rocket  bodies,  and  junk  is  disproven  by  costly  collisions.  For  example,  the  col¬ 
lision  of  the  Iridium  33  satellite  with  the  defunct  Russian  Cosmos  2251  satellite  in 
February  2009  was  not  predicted  ahead  of  time.  The  collision  was  only  known  to 
have  occurred  after  Iridium  33  become  unresponsive  and  US  space  operators  began 
tracking  a  debris  cloud  which  originated  at  Iridium  33 ’s  previous  position.  While  the 
two  satellites  were  flagged  by  the  Center  for  Space  Standards  and  Innovation  (CSSI) 
as  having  a  close  approach  of  just  over  half  of  a  kilometer,  there  were  other  closer 
approaches  predicted  that  day  (4). 

The  problem  of  identifying  probable  collisions  with  sufficient  lead-time  and  ac¬ 
curacy  cannot  be  solved  with  today’s  methods.  Enter  the  KAM  Theorem.  First  pro¬ 
posed  by  Andrey  Kolmogorov  in  1954,  and  later  extended  and  proved  by  Vladimir 
Arnold  in  1963  and  Jurgen  Moser  in  1962,  the  KAM  Theorem  potentially  offers  a 
new  method  of  attacking  this  problem.  The  KAM  Theorem,  when  applied  to  Earth 
orbiting  satellites,  essentially  transforms  the  physical  variable  state  description  of  the 
system  to  KAM  Torus  variables  in  which  three  coordinates  increment  linearly  with 
time  and  three  momenta  remain  constant.  The  KAM  Torus  can  be  visualized  as  a 
multi-dimensional  donut  on  which  the  satellite  traverses  the  surface.  The  KAM  Theo¬ 
rem  may  also  provide  a  replacement  to  the  current-day  two- line  element  sets  now  used 
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to  publish  a  satellite’s  ephemerides  which  are  subsequently  used  in  many  applications 
including  satellite  collision  predictions.  The  replacement  is  only  now  possible  due  to 
today’s  computing  power  combined  with  the  new  application  of  the  KAM  Theorem 
and  its  ability  to  consider  the  entire  Earth’s  geopotential  not  as  a  perturbation,  but 
as  a  fundamental  component  of  the  orbit  description. 

The  study  of  orbital  perturbations  follows  one  of  two  approaches:  general  pertur¬ 
bations  and  special  perturbations.  General  perturbation  theory  is  mainly  analytical 
and  seeks  to  understand  how  on  average  various  perturbations  affect  a  satellite’s  or¬ 
bital  elements.  General  perturbation  theory  will  offer  understanding  to  the  orbital 
mechanist,  but  not  the  precision  knowledge  required  for  today’s  applications  or  for 
predicting  and  preventing  satellite  collisions.  Special  perturbation  theory  well  reflects 
its  name  -  it  involves  numerically  integrating  the  equations  of  motion  to  determine 
the  specific  solution  of  a  certain  satellite  and  its  initial  conditions.  Change  the  ini¬ 
tial  conditions  and  the  numerical  integration  will  have  to  be  re-performed  for  a  new, 
unique  solution.  Considering  the  incredible  number  of  artificial  satellites  orbiting  the 
Earth,  the  special  perturbations  approach  to  predicting  and  preventing  satellite  colli¬ 
sions  is  inefficient  if  not  impossible.  In  the  search  for  new  methods  KAM  Theory  may 
provide  the  answer. 

As  more  nations  around  the  world  gain  homegrown  access  to  space,  the  num¬ 
ber  of  artificial  satellites  purposefully  put  into  space  to  orbit  the  Earth  will  increase 
dramatically.  But  along  with  the  purposeful  increase  will  come  the  passive  increase. 
As  the  satellite  count  increases,  so  too  does  the  probability  for  collisions.  As  an  illus¬ 
tration,  ponder  this  message  from  T.S.  Kelso,  posted  on  Twitter  on  January  4,  2011 
(10):  UTLE  data  released  for  76  more  pieces  of  Cosmos  2251  debris,  bringing  total  to 
1,423  with  only  75  pieces  having  decayed .”  In  less  than  two  years,  1,348  known  pieces 
of  debris  from  the  Iridium  and  Cosmos  collision  remain.  Not  considering  the  debris 
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yet  to  be  identified,  these  1,348  pieces  must  be  tracked;  however,  just  because  debris  is 
tracked  does  not  mean  it  cannot  or  will  not  collide  with  other  satellites  in  the  future. 

1.2  Approach 

In  this  work,  the  author  has  attempted  to  quantify  how  air  drag  and  third-body 
lunar  perturbations  affect  orbital  motion  near  a  reference  KAM  Torus.  Expressions 
for  the  average  effects  to  the  torus  coordinates  and  momenta  were  determined.  The 
satellite’s  perturbed  state  was  numerically  integrated  and  transformed  into  torus  co¬ 
ordinates  and  momenta  for  comparison  to  the  nearby  reference  torus.  A  simple  at¬ 
mosphere  model  was  utilized  for  drag  calculations.  The  perturbing  accelerations  due 
to  air  drag  and  lunar  third-body  were  incorporated  in  the  numerical  integration.  For 
simplicity,  the  moon  was  assumed  to  orbit  in  the  equatorial  plane. 

1.3  Problem  Statement 

Since  the  KAM  Theorem  is  based  on  a  nearly  integrable  Hamiltonian  system, 
it  is  necessary  to  show  that  application  of  the  theorem  to  low-earth  orbiting  satellites 
is  not  precluded  by  air  drag,  which  is  a  non-conservative  force,  or  lunar-third  body 
effects.  The  deviation  of  the  perturbed  motion  from  the  reference  torus  should  be 
quantified  and  understood.  If  the  deviation  is  non-negligible,  a  method  should  be 
proposed  for  accounting  for  these  perturbing  effects. 
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II.  Background 


The  simple  case  of  an  object  orbiting  a  larger  body  is  described  as  a  central 
force  problem  in  which  the  force  is  due  to  gravity.  For  a  central  force  problem 


F  =  F,r-  (1) 

where  Fg  is  defined  in  Equation  2  and  f,  is  the  unit  position  vector: 

GMm 
* a  =  ~2 

where  G  is  the  gravitational  constant,  M  is  the  mass  of  the  larger  body,  and  m  is  the 
mass  of  the  smaller  body. 

This  version  of  the  central-force  problem  is  often  referred  to  as  the  “Keplerian” 
problem  and  results  in  the  familiar  orbital  solutions  in  which  the  state  of  an  orbiting 
mass  is  known  for  all  time.  Especially  in  satellites  orbiting  the  Earth,  it  can  be  readily 
seen  that  reality  quickly  diverges  from  the  Keplerian  solution. 

2.1  Perturbations  to  the  Keplerian  Problem 

Over  the  years,  orbital  mechanists  have  shown  that  many  forces  perturb  the 
solutions  predicted  by  the  two-body  problem.  Earth  orbiting  satellites  are  subject  to 
the  gravitational  potential  differences  of  the  non-spherical  Earth.  Spacecraft  orbiting 
the  Earth  at  low  altitudes  are  susceptible  to  force  resulting  from  air  drag.  Third- 
body  effects  from  the  Moon  also  perturb  the  Keplerian  solution.  Other  sources  of 
perturbation  include  the  Sun’s  gravity  and  solar  radiation  pressure,  but  these  will  not 
be  discussed  in  this  work. 
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2.1.1  Geopotential.  The  assumption  of  a  “point-mass”  is  made  in  the  deriva¬ 
tion  of  Keplerian  motion.  For  accurate  applications,  the  point-mass  assumption  cannot 
be  made  as  the  Earth  is  not  perfectly  spherical.  The  shape  of  the  Earth  is  actually 
that  of  a  spheroid  -  the  shape  that  is  obtained  by  rotating  an  ellipse  about  one 
of  its  axes.  This  mass  distribution  results  in  a  gravity  potential  function  that  can  be 
described  using  zonal,  sectorial,  and  tesseral  harmonics.  Zonal  harmonics  are  indepen¬ 
dent  of  longitude  and  include  the  famous  oblateness  term,  “  J2.”  Sectorial  harmonics 
are  independent  of  latitude  and  resemble  orange  slices.  Finally,  tesseral  harmonics  are 
dependent  on  both  longitude  and  latitude  and  resemble  a  chessboard  placed  over  the 
surface  of  the  Earth  (7). 

The  Earth’s  geopotential  function  is  given  by  Wiesel  (23)  as 


oo  n  /  \  —n 

V  (r,  6,<j>)  =  ~  '  f  )  p™  (cos  9)  (Cnm  cos  m(f>  +  Snm  sin  m<j>)  (3) 

n= 0  m= 0  '  e  ' 

where  r,  9,  and  <p  are  spherical  coordinates;  R,  is  the  radius  of  the  Earth;  //  is  the 
Earths  gravitational  parameter;  Cnm  and  Snrn  are  constants  specifying  the  shape  of  the 
gravitational  field;  and  P™  are  associated  Legendre  polynomials.  When  m  =  n  =  0, 
the  geopotential  function  degenerates  into  the  Keplerian  potential  given  by 

V(r)  =  -t  (4) 

r 

Ideally,  one  could  use  the  entire  geopotential  to  derive  the  perturbations  to  the 
Keplerian  motion  of  a  satellite.  Unfortunately,  developing  the  perturbations  analyt¬ 
ically  for  just  the  second  zonal  harmonic  (to  =  0,  n  =  2)  is  very  complex.  On  the 
other  hand,  since  the  second  zonal  harmonic  is  by  far  the  largest  contributor  to  non- 
Keplerian  effects,  these  perturbations  to  the  Keplerian  solution  can  be  isolated  and 
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well  understood.  The  average  rate  of  change  of  the  orbital  elements  due  to  the  second 
zonal  harmonic  will  not  be  derived  here  but  are  listed  by  Schaub  (16): 


d^  =  \J2ni^)  (3 cos2  f  —  l)  (10) 

where  a  is  the  semi-major  axis,  e  is  the  eccentricity,  i  is  the  inclination,  fl  is  the  right 
ascension  of  the  ascending  node  (RAAN),  J>  is  a  constant,  n  is  the  mean  motion,  Re  is 
the  radius  of  the  Earth,  p  is  the  semi-latus  rectum,  c 0  is  the  argument  of  perigee  (AP), 
and  M  is  the  mean  anomaly.  Note  that  Equations  5,  6,  7,  8,  9,  and  10  are  the  secular 
rates  and  do  not  exhibit  short  or  long  period  oscillations. 

2.1.2  Air  Drag  Perturbations.  The  Earth’s  atmosphere  extends  to  approxi¬ 
mately  1500  km  altitude.  Satellites  orbiting  in  the  regime  below  1000  km  altitude  are 
subject  to  a  non- negligible  force  due  to  air  drag  that  is  a  function  of  cross-sectional 
area,  mass,  drag  coefficient,  velocity,  and  air  density  (6).  Because  drag  is  a  non¬ 
conservative  force,  the  satellite  will  continuously  lose  energy  when  exposed  to  the 
drag  force.  Elliptical  orbits  will  circularize  as  the  apogee  altitude  is  continuously  re¬ 
duced  due  to  the  energy  loss.  Once  the  orbit  has  been  circularized,  the  net  effect  is  a 
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reduction  of  orbital  altitude  until  the  satellite’s  velocity  is  reduced  to  the  point  that 
it  cannot  stay  in  orbit. 

2. 1.2.1  Effects  of  Air  Drag  on  Orbital  Elements.  Schaub  describes  (16) 
inserting  the  acceleration  due  to  air  drag  into  Gauss’  variational  equations  to  arrive  at 
average  rates  of  change  to  the  classical  orbital  elements.  Gauss  variational  equations 
are  listed  below: 


de 

dt 


da 

dt 


2  a2V 
d 


ay 


^  sin  v  an  +  2  (e  +  cos  v)  ay^j 


di  r  cos  (lo  +  v) 

—  =  — .  —  ah 

dt  ^/pa  (1  —  e2) 

dfl  r  sin  (uj  +  v) 
dt  sin  i^fpfa  (1  —  e2) 


(11) 

(12) 

(13) 

(14) 


doj  If  (  r\  \  r  sin  (lu  +  v)  cos  i 

—  =  —  (^-an  y2e  +  -J  cos  v  +  2ay  sin  uj - y  ^  ah 


\f  da  (1  -  e2) 


dM 


=  n  + 


VT ~ 


—a„  cos  v  —  2(1 


2 

e r 


ay  sm  v 


(15) 


(16) 


dt  '  eV  \a  \  a  (1  —  e2) 

where  V  is  the  velocity;  ay,  ah,  and  an  are  the  accelerations  in  the  local- vertical-local- 
horizontal  (LVLH)  coordinate  frame;  ft  is  the  Earth’s  gravitational  parameter;  r  is  the 
radius;  and  v  is  the  true  anomaly. 


Air  drag  only  occurs  in  the  velocity  direction,  ah  =  0  and  an  =  0,  leaving  only 
ay  =  af)  =  —V2/2B*p  which  when  placed  in  the  variational  equations  yields  (16): 


da 

dt 


de 

dt 


—B*pV  (e  +  cos  v) 


(17) 
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(18) 


(19) 


di 

dd=° 

dn  =  0 

dt 

A  =  -BV 

dt  e 


dM 


=  n  +  B*  pV 


a/1  -  e2 


l  +  eJ 


sm  v 


(20) 

(21) 

(22) 


dt  '  r  e  ^  '  a  (1  —  e2) 

where  B*  is  the  ballistic  drag  coefficient  and  p  is  the  atmospheric  density  at  altitude. 
These  equations  give  us  insight  into  the  effect  of  drag  on  a  satellites  orbit.  Quite 
obviously,  the  semi-major  axis  will  shrink  but  the  other  effects  really  require  the 
parameters  to  be  evaluated  as  an  orbit  is  propagated  to  gain  understanding. 


Note  that  as  eccentricity  approaches  a  zero  value,  these  equations  are  invalid 
since  the  argument  of  perigee  and  right  ascension  of  the  ascending  note  are  undefined. 
The  eccentricity  in  the  denominator  of  two  of  the  equations  would  present  a  problem 
as  well. 


2. 1.2. 2  Atmosphere  Model.  The  properties  of  the  atmosphere  are  di¬ 
verse  by  location  and  time.  A  simplistic  model  that  only  considers  variation  by  altitude 
is  the  exponential  model  in  which  the  density  of  the  atmosphere  is  simply  a  function 
of  a  reference  density  and  height: 


—  (r  —  rn) 

p(r)  =  poe~^  (23) 

where  po  is  the  reference  density  at  some  reference  position,  ro,  and  H  is  the  scale 
height. 

For  the  purposes  of  this  research,  a  more  realistic  model  must  be  considered  be¬ 
cause  above  86  km  altitude,  the  atmosphere  is  highly  dynamic  (15).  The  model  chosen 
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for  this  research  is  detailed  by  Regan  and  Anandakrishnan  (R&A).  The  atmosphere 
is  considered  to  be  broken  into  layers  of  atmospheric  strata.  The  boundary  between 
each  stratum  is  considered  to  occur  when  the  atmosphere  is  isothermal.  R&A  provide 
a  TRUEBASIC  algorithm  that  calculates  the  atmospheric  properties  of  temperature, 
pressure,  and  density,  for  a  given  altitude.  The  R&A  Standard  Atmosphere  is  a  hybrid 
model  in  which  the  U.S.  1976  Standard  Atmosphere  is  used  for  altitudes  from  0  to  86 
km  and  the  U.S.  1962  Standard  Atmosphere  is  used  for  altitudes  above  86  km.  R&A 
created  this  hybrid  model  for  sake  of  algorithm  simplicity.  Above  86  km,  the  1976 
model  contains  two  strata  in  which  the  lapse  rates  are  not  linear  but  rather  elliptical 
and  exponential.  R&A  make  the  argument  that  while  these  rates  could  be  calculated, 
their  minor  difference  from  the  1962  model  are  easily  overcome  by  error  introduced 
by  incoming  solar  radiation. 

The  author  has  translated  the  R&A  Standard  Atmosphere  TRUEBASIC  algo¬ 
rithm  into  C++  for  compatibility  with  Wiesel’s  software  (21).  The  translated  algo¬ 
rithm  provides  only  density  as  it  is  the  only  output  required  for  this  research.  Figure 
1  shows  the  density  of  the  atmosphere  from  0  to  700  km  per  the  adapted  R&A  algo¬ 
rithm. 

The  algorithm  was  validated  by  comparing  its  output  to  the  altitude/density 
profile  of  the  U.S.  1962  &  1976  Standard  Atmospheres  published  by  National  Aero¬ 
nautics  and  Astronautics  Administration  (NASA).  Figure  2  shows  the  1962  &  1976 
U.S.  Standard  Atmosphere  and  three  earlier  U.S.  Standard  and  model  atmospheres. 

Even  more  advanced  models  must  consider  many  additional  factors  including: 
longitude  and  latitude,  time  of  day,  time  of  the  year,  solar  cycle,  and  solar  storms.  For 
the  purposes  of  this  research,  models  beyond  the  complexity  of  the  chosen  model  will 
not  be  utilized  since  the  intent  is  to  understand  how  air  drag  modifies  the  KAM  Torus 
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Figure  1:  Atmospheric  Density  Plot  Generated  From  Translated  R&A  Standard  At¬ 
mosphere  Algorithm 

coordinates  and  momenta  in  general,  rather  than  provide  precise  state  predictions  in 
a  specific  scenario. 


2.1.3  Moon  Perturbations.  The  third-body  problem  has  been  treated  by 
numerous  distinguished  orbital  mechanists  (3),  such  as  L.  Euler,  C.  E.  DeLaunay,  H. 
Poincare,  G.  W.  Hill,  and  E.  W.  Brown  as  well  as  D.  Brouwer  (3),  D.  Vallado  (17), 
Schaub  (16),  Chao  (6)  and  many  more.  Due  to  its  difficulty,  the  most  insightful  results 
have  been  garnered  when  the  three  bodies  are  restricted  to  in  plane  motion. 

Chao  (6)  writes  that  the  disturbing  function  for  third-body  motion  is  given  as 


r3 


r  cos  S 

r3 


(24) 


11 


1000 

900 

800 

700 

—  600 
E 
jc 

-e 
? 

x  500 
u 

• 

300 

200 

100 

10"1S  10',s  10'"  10  s  10'7  10  s  lO'*  10'1  0  101 

Figure  2:  1962  &  1976  U.S.  Standard  Atmosphere  (14) 

where  /i3  is  the  gravitational  constant  of  the  third-body,  r3  is  the  position  vector 
magnitude  of  the  third-body,  r  is  the  position  vector  magnitude  of  the  satellite,  and 
S  is  the  angle  between  the  two  position  vectors  r  and  r3. 

Kaufman  (9)  showed  that  when  r/r3  <<  1,  the  square-root  term  in  Equation 
24  can  be  expanded  with  higher  order  terms  discarded.  Kaufman  then  applied  the 
method  of  averaging  to  remove  the  short  period  terms,  resulting  in  a  new  expression 
for  the  disturbing  function: 
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where  A  =  P  •  e3  and  B  =  Q  •  e3;  a  is  the  semi-major  axis  and  e  is  the  eccentricity 
of  the  satellite;  a3  is  the  semi-major  axis  of  the  third-body,  r3  is  the  position  vector 
magnitude  of  the  third-body,  and  n3  is  the  mean  motion  of  the  third  body;  and 


cos  D  cos  u j  —  sin  Q  sin  u  cos  i 
sin  Q  cos  to  +  cos  sin  u j  cos  i 
sin  uj  sin  i 

—  cos  sin  u  —  sin  D  cos  uo  cos  i 
Q  =  —  sin  fl  sin  u  +  cos  fl  cos  u  cos  i 

cos  uj  sin  i 

cos  fl3  cos  u3  —  sin  fl3  sin  u3  cos  i3 
e3  =  sin  fl3  cos  u3  +  cos  fl3  sin  u3  cos  i3  (28) 

sin  u3  sin  i3 

where  e3  is  the  unit  postion  vector  of  the  third-body,  i3  is  the  inclination  of  third- 
body,  u3  is  the  argument  of  latitude  of  the  third-body,  and  fi3  is  the  RAAN  of  the 
third-body. 

Chao  (5)  then  expanded  the  disturbing  function  per  the  classical  orbital  ele¬ 
ments.  The  final  results  (6)  are  equations  for  the  rate  of  change  of  the  classical  orbital 
elements  of  the  satellite  that  are  being  perturbed  by  the  third-body.  The  equations 
are  too  lengthy  to  cover  here  but  they  are  heavily  dependent  on  the  AP  and  RAAN 
of  both  the  third-body  and  the  satellite.  Due  to  assumptions  that  will  be  made  to 
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simplify  the  Moon  model  for  this  research  -  chiefly  that  the  Moon’s  orbit  is  circular 
about  the  Earths  equator  (AP  and  RAAN  are  undefined)  -  these  equations  will  not 
provide  results  for  comparison  in  this  research. 

Vallado  (17)  summarizes  the  effects  on  a  satellite  due  to  third-body  perturba¬ 
tions  as  periodic  for  semi-major  axis,  eccentricity,  and  inclination  and  secular  and 
periodic  for  mean  anomaly,  RAAN,  and  AP.  An  approach  utilizing  KAM  Theory  will 
be  discussed  in  §3.3  -  the  results  will  be  compared  against  this  summary. 

2.2  KAM  Theory 

The  KAM  Theorem  was  developed  and  proved  over  the  period  of  about  a  decade 
by  Andrey  Kolmogorov  (11),  Vladimir  Arnold  (1),  and  Jurgen  Moser  (13).  Essentially, 
the  KAM  Theorem  states  that  the  solutions  of  a  slightly  perturbed  Hamiltonian 
system  will  exist  on  the  surface  of  a  torus. 

KAM  theory  is  born  out  of  a  special  subclass  of  Hamiltonian  dynamics  that 
exhibit  periodic  motion  at  certain  fundamental  frequencies.  The  consequence  of  being 
characterized  as  a  periodic  Hamiltonian  system  is  that  the  physical  variables  can 
be  transformed  into  coordinate-momentum  pairs  in  which  the  momenta  are  constant 
and  the  coordinates  are  angles  that  increment  linearly  with  time.  This  becomes  one 
of  the  primary  benefits  of  describing  an  orbital  system  as  a  KAM  Torus.  In  current 
application  of  the  KAM  theory  to  Earth  orbiting  systems,  the  transformation  results 
in  three  system  coordinates  that  increment  linearly  with  time 

Q  =  Qo  T  kV  (^9) 

where  Q  is  the  coordinate  vector,  Q0  is  the  initial  values  of  the  coordinate  vector, 
and  u>  is  the  basis  frequency  vector. 
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In  addition  to  the  coordinates,  three  constant  momenta,  P,  result  from  the 
transformation  to  the  KAM  Torus.  Predicting  the  future  state  of  the  system  becomes 
trivial  in  this  scenario.  The  application  of  KAM  Theory  to  Earth  orbiting  systems 
offers  another  tempting  benefit  to  the  orbital  mechanic:  the  Earth’s  geopotential  is 
incorporated  in  the  KAM  Torus  to  an  arbitrary  order,  thus  converting  geopotential 
perturbations  in  physical  variables  into  standard  torus  behavior. 

The  application  of  KAM  Theory  to  Earth  satellites  had  remained  relatively  un¬ 
touched  until  Wiesel  first  argued  that  at  least  some  Earth  orbits  can  be  described  by 
KAM  Tori  (20).  The  evidence  confirming  this  argument  was  first,  that  the  Earth  satel¬ 
lite  torus  contained  only  three  fundamental  frequencies  as  predicted  by  the  Hamilton- 
Jacobi  theorem;  and  second,  that  the  motion  of  the  satellite  could  be  described  via  a 
Fourier  series  expansion.  Wiesel  also  showed  that  the  two-body  problem  in  the  earth- 
centered-earth- fixed  (ECEF)  frame  and  the  simplified  general  perturbations  4  (SGP4) 
model  can  be  considered  to  be  torus  models  as  well  (22).  In  the  same  paper,  Wiesel 
also  showed  that  the  KAM  Torus  is  an  effective  perturbation  theory  for  an  Earth 
satellite  since  the  entire  Earth  geopotential  is  incorporated  in  the  torus. 

Several  doctoral  and  graduate  students  of  Wiesel  have  also  contributed  to  the 
body  of  Earth  satellite  KAM  Tori  research.  B.  Little  attempted  to  fit  tori  to  actual 
data  from  the  Jason- 1  and  Gravity  Recovery  and  Climate  Experiment  (GRACE) 
satellites  (12).  He  determined  that  Jason-1  does  indeed  lie  on  a  torus.  On  the  other 
hand,  he  was  unable  to  fit  GRACE  to  a  torus  and  proposed  that  perturbation  due  to 
air  drag  may  have  been  responsible.  C.  Craft  investigated  the  applicability  of  KAM 
Theory  to  Earth  orbiting  satellite  formations  (8) .  He  determined  that  formations  with 
large  physical  separations  have  proportionally  large  secular  drift  which  may  preclude 
the  use  of  KAM  Theory  for  formation  flight  applications.  Additionally,  he  showed 
that  smaller  formations  in  which  the  satellites  have  small  separations  in  the  torus 
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coordinates  result  in  drift  rates  easily  countered  by  current-day  technology  such  as 
electric  propulsion.  R.  Bordner  fit  tori  to  Global  Positioning  System  (GPS)  satellite 
orbital  data  and  found  that  the  resonant  nature  of  the  GPS  orbit  resulted  in  the 
smallest  system  frequency  being  masked  by  the  fact  that  an  entire  period  did  not 
occur  within  the  orbital  data  set  (2).  This  trouble  led  Bordner  to  develop  a  new 
procedure  of  fitting  the  frequency  clusters  within  the  orbital  data  to  the  analytical 
form  of  the  truncated,  continuous,  Fourier  transform  (ATCFT)  to  form  the  KAM 
Torus. 
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III.  Methodology 


The  overall  approach  to  understanding  how  motion  near  a  KAM  Torus  is  per¬ 
turbed  by  air  drag  and  the  moon  will  be  to  propagate  the  perturbed  orbit  via  nu¬ 
merical  integration,  propagate  the  non-perturbed  torus  orbit  by  updating  the  torus 
coordinates  as  described  in  §2.2,  and  then  compare  the  two  resulting  states  in  terms 
of  torus  coordinates  and  momenta.  These  state  differentials  will  then  be  transformed 
back  to  physical  coordinates  to  aid  in  understanding  the  type  and  quantity  of  error. 

Specifically,  the  remainder  of  the  research  will  consider  and  compare  three  cases: 

1.  Baseline :  comparison  of  unperturbed  numerically  integrated  motion  to  reference 
torus  motion.  Any  error  between  the  two  should  be  due  to  linearizations  and 
numerical  truncation 

2.  Drag:  comparison  of  drag  perturbed  numerically  integrated  motion  to  reference 
torus  motion 

3.  Lunar:  comparison  of  third-body  lunar  perturbed  numerically  integration  mo¬ 
tion  to  reference  torus  motion 

Discussion  and  figures  in  the  remaining  sections  of  this  work  will  refer  to  these 
cases  as  “Three  Cases”  or  “Baseline,  Drag,  and  Lunar  Cases.”  Wiesel  has  developed 
a  C++  software  package  (21),  that  evaluates  the  motion  of  a  satellite  near  a  reference 
torus.  The  inner  workings  of  the  software  will  be  described  in  §§3.2. 1-2.  The  software 
was  specifically  designed  to  examine  the  motion  of  an  Earth  satellite  with  initial 
coordinates  and/or  momenta  offsets  from  the  reference  torus.  The  author  has  modified 
this  software  to  evaluate  perturbed  motion  near  a  reference  torus  in  the  following  ways: 

•  Addition  of  standard  atmosphere  object  for  calculating  atmospheric  density  at 
altitudes  0-700  km 
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•  Calculation  of  position  and  velocity  error  for  total  or  single  state  differences 
between  perturbed  numerically  integrated  motion  and  reference  torus  motion 

•  Calculation  of  rates  of  change  of  torus  coordinates  and  momenta 

•  Option  to  execute  software  in  three  additional  configurations: 

—  Motion  perturbed  by  air  drag 
—  Motion  perturbed  by  lunar  third-body  effects 
—  Motion  perturbed  by  both  air  drag  and  lunar-third  body 

A  specific  description  of  how  the  software  was  modified  to  incorporate  perturb¬ 
ing  accelerations  due  to  air  drag  and  the  Moon’s  third-body  effects  is  presented  in 
§§3.2. 2. 2-3. 

The  author  performed  data  analysis  and  plot  generation  using  the  open  source 
scripting  language  Python  extended  by  the  following  open  source  packages: 

•  NumPy  -  numerical  package  that  offers  n-dimensional  array  object  as  well  as 
many  numerical  functions  for  operations  on  arrays 

•  SciPy  -  scientific  package  that  offers  functions  for  integration,  optimization, 
linear  algebra,  fourier  transforms,  etc 

•  Matplotlib  -  plotting  package  with  functionality  similar  to  common  numerical 
plotting  tools 

3.1  Obtaining  the  Reference  KAM  Torus 

3.1.1  Method  for  Finding  the  KAM  Torus.  While  several  methods  exist  for 
finding  the  KAM  Torus  of  an  Earth  orbiting  satellite,  the  following  method,  estab¬ 
lished  by  Wiesel  and  Bordner  (2)  was  used.  First,  an  orbit  with  an  initial  state  is 
chosen  and  numerically  integrated,  including  only  geopotenial  perturbations,  forward 
and  backward  in  time  for  six  months  in  both  directions.  Next,  the  raw  position  data 
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from  the  numerical  integration  of  the  orbit  is  transformed  into  a  spectral  plot.  The 
spectral  plot  is  then  sectioned  into  groups  of  frequency  lines  where  the  analytical  form 
of  the  truncated,  continuous,  Fourier  transform  (ATCFT)  is  used  to  fit  the  frequency 
lines  and  determine  the  Fourier  coefficients  (2).  Bordner  terms  this  as  the  “Frequency 
Cluster-Based  Approach.” 

3.1.2  KAM  Tori  Utilized  in  Research.  The  author  utilized  two  tori  in  this 
research.  Torus  #l’s  classical  orbital  element  characterization  is  listed  in  Table  1. 
The  classical  orbital  elements  were  converted  to  the  associated  position  and  velocity 
vectors  and  integrated  forward  and  backward  in  time  by  6  months  for  one  year  of 
total  data.  The  torus  was  then  created  using  the  approach  developed  by  Bordner  (2) 
as  described  in  §3.1.1. 

The  second  torus  was  provided  by  Yates  (24).  He  obtained  high  precision  orbital 
data  of  the  International  Space  Station  (ISS).  He  then  passed  the  data  through  a 
Bayes  filter  several  times  in  order  to  condition  the  data  such  that  it  was  conducive 
to  matching  with  a  KAM  Torus.  The  torus  was  finally  constructed  as  described  in 
§3.1.1  with  the  exception  that  the  conditioned  ISS  orbital  data  was  used  in  place  of 
the  numerical  integration  data.  The  classical  orbital  element  characterization  of  the 
ISS  orbit  is  listed  in  Table  1: 

Table  1:  Classical  Orbital  Elements  of  Orbits  Represented  by  Tori 


Torus 

a  [DU] 

e 

i° 

n° 

w° 

M° 

#1 

1.1 

0.05 

30.00 

261.7 

141.4 

88.4 

#2  (ISS) 

1.05 

0.001 

51.6 

300.2 

112.7 

322.4 

Though  two  tori  were  available  for  analysis,  the  majority  of  this  work  consid¬ 
ers  only  Torus  #1  because  the  results  for  each  torus  are  very  similar.  Results  are 
assumed  to  refer  to  Torus  #1  unless  indicated  otherwise.  Note  that  both  tori  have 
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semi-major  axes  that  are  within  the  Earth’s  atmosphere  and  will  be  affected  by  air 
drag  (important  if  the  effects  of  air  drag  on  tori  are  to  be  analyzed). 

3.2  Numerical  Integration  of  Perturbed  Motion  Near  the  Reference 
Torus 

3.2.1  Extracting  Initial  State  from  the  Torus.  The  initial  conditions  in  phys¬ 
ical  coordinates  must  be  established  from  the  torus  via  multiple  Fourier  series  (22): 

q  [Cj  cos  (j  •  Q)  +  Sj  sin  (j  •  Q)]  (30) 

j 

where  Cj  and  Sj  are  vector  coefficients  with  a  vector  summation  index  j,  and  Q 
are  the  torus  coordinates.  Next,  the  initial  momenta  must  be  determined  through 
manipulation  of  the  Lagrangian: 

dC 

where  p  are  the  physical  momenta,  C  is  the  Lagrangian,  and  q  are  the  inertial  time 
derivatives  of  the  physical  coordinates. 

3.2.2  Equations  of  Motion.  The  equations  of  motion  will  be  addressed  in 
three  parts.  The  primary  effort  considers  the  two-body  problem  perturbed  by  the 
geopotential.  Once  this  effort  is  properly  posed,  it  will  be  extended  to  air  drag  and 
Moon  third-body  motion  perturbations. 

3.2.2. 1  Keplerian  Motion  Plus  Geopotential.  In  the  Cartesian  space 
rotating  with  the  Earth,  the  ECEF  coordinate  system,  the  Hamiltonian  for  an  orbiting 
satellite  is  derived  as  any  other  Hamiltonian,  beginning  with  Equation  32. 
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(32) 


H  =  -  C 

i 

where  Pi  are  momenta,  qt  are  the  coordinate  velocities,  C  is  the  Lagrangian,  and  H  is 
the  Hamiltonian. 

The  Lagrangian  is  calculated  by  subtracting  a  system’s  potential  energy  from  a 
term  representing  the  system’s  kinetic  energy  as  represented  by  Equation  33. 


C  =  T  -V  (33) 

where  T  is  the  kinetic  energy  term  and  V  is  the  potential  energy  term. 

The  energy  term  can  easily  be  calculated  via  Equation  36.  Start  with  the  satel¬ 
lite’s  position  expressed  in  coordinates  x,  y,  and  z  in  the  ECEF  frame.  The  inertial 
derivative  of  the  position  can  be  obtained  with  the  help  of  the  Transport  Theorem 
and  is  shown  by  Equation  34: 


r 


EC  I 


r  =i 


ECEF 


r  +  x  r 


(34) 


which  when  expanded  yields 


x  -  yu 

Qi 

y  +  xuj ® 

— 

Q2 

i 

Q3 

(35) 


Next,  now  that  the  velocity  is  known  in  terms  of  the  coordinates,  the  kinetic 
energy  per  unit  mass  can  be  calculated: 


T  =  ^V2  =  ^[(x-  yu®)2  +  (y  +  xw©)2  +  z2] 


(36) 
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Now  the  Lagrangian  becomes 


C  =  7.  [(£  -  ! Ms)  +  (</  +  +  22]  -  V(x,  a,  2) 


(37) 


Next,  to  calculate  the  momenta,  Equation  (34)  is  utilized: 


dC 


The  momenta  are 


<9q 


P 


x  -  yco © 
V  +  xujQ 


(38) 


(39) 


For  use  in  calculation  of  the  Hamiltonian,  Equation  39  is  rearranged  to  solve 
for  the  rate  of  change  of  the  coordinates: 


Pi  -  yu® 

q  =  p2  +  xuje 

P3 

where  pi,  p2,  and  p3  are  the  physical  momenta. 


(40) 


Now  utilizing  Equations  32,  37,  39,  and  40  the  Hamiltonian  is  found  to  be 


'H  =  \  (pi  +  p\  +  p\)  +  ^®ypi  -  u®xp2  +  V  (41) 

where  H  is  the  Hamiltonian,  x,  y,  and  z  are  the  coordinates,  ui®  is  the  Earth’s  angular 
velocity,  and  V  is  the  Earth’s  geopotential  as  described  in  Equation  3. 
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The  coordinates  and  momenta  are  the  system  state: 


X  = 


Qi 

X 

Q2 

y 

*73 

z 

Pi 

Vx 

P2 

Vy 

P3 

V 

(42) 


where  q\,  q-2,  and  q?,  are  the  physical  coordinates  and  Vx,  Vy,  and  Vz  are  the  physical 
coordinate  velocities. 


Since  the  ECEF  coordinate  system  is  being  used,  the  Earth’s  gravitational  po¬ 
tential,  V,  is  not  a  function  of  time  but  of  only  the  satellite’s  position.  Therefore, 
Hamilton’s  equations  can  easily  be  used  to  determine  the  time  rate  of  change  of  the 
coordinates  and  momenta: 


dH 


Pi  = 


dH 

dqi 


(43) 

(44) 


Because  the  coordinates  are  the  position  and  the  momenta  are  the  velocity, 
Equations  43  and  44  can  be  used  to  determine  how  position  and  velocity  of  the  Earth 
satellite  are  changing  with  time: 


_  dH 

~  dVx 


Vx  + 


(45) 


y  = 


dH 

Wy 


Vy  ~ 


(46) 
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(47) 


z  = 


dU 

dVz 


K 


V, 


m  _  ov 


(48) 


dV 

dy 


(49) 


m  _  _ov 

dz  dz 


(50) 


Note  that  while  these  equations  of  motion  were  derived  with  the  Hamiltonian 
(i.e.  conservative  system),  accelerations  due  to  non-conservative  forces  or  additional 
conservative  forces  can  be  added  to  Vx,  Vy,  and  14  with  ease.  This  approach  will  be 
taken  in  §§3.2. 2. 2-3  when  specifically  addressing  air  drag  and  the  third-body  Moon 
problem. 

Now  that  the  equations  of  motion  are  established,  numerical  integration  can  be 
performed.  The  Hamming  fourth  order  predictor /corrector  numerical  integrator  was 
utilized.  While  the  equations  of  motion  only  contain  accelerations  due  to  conservative 
forces,  the  accuracy  of  the  numerical  integration  can  be  checked  by  evaluating  the 
Hamiltonian  at  each  time  step.  Once  air  drag  or  lunar  effects  are  added  this  check 
will  not  be  valid. 


3. 2. 2. 2  Air  Drag.  The  incorporation  of  air  drag  in  the  equations  of 
motion  is  as  simple  as  modifying  Equations  48,  49,  and  50  to  include  the  following 
vector  term: 
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(51) 


V  =  -\pv^B‘  =  -\pWBT 

where  p  is  the  density  of  the  atmosphere,  V  is  the  velocity  of  the  satellite  relative  to 
the  atmosphere,  V  is  the  velocity  vector  in  the  ECEF  frame,  and  B*  is  the  ballistic 
coefficient. 

The  ballistic  coefficient  is  dependent  on  the  drag  coefficient,  cross-sectional  area 
in  the  velocity  direction,  as  well  as  the  mass.  Note  that  the  dimensions  of  the  ballistic 
coefficient  must  be  area  per  mass  for  Equation  51  to  be  valid.  Many  Earth  orbiting 
satellites  have  rotational  motion  about  their  center  of  mass,  which  results  in  the 
values  of  the  drag  coefficient  and  cross-sectional  area  varying  with  time.  For  the 
purposes  of  this  research,  the  ballistic  coefficient  will  be  assumed  to  be  constant. 
Table  2  summarizes  ballistic  coefficients  for  several  satellites: 

Table  2:  Ballistic  Coefficients  of  Various  Satellites  (18) 


Satellite 

Shape 

Max. 

Ballistic 

Coef. 

[kg/m2] 

Min. 

Ballistic 

Coef. 

[kg/m2] 

Max. 

Ballistic 

Coef. 

[m2/kg] 

Min. 

Ballastic 

Coef. 

[m2/kg] 

Oscar- 1 

Box 

42.8 

16.7 

0.0234 

0.0599 

Intercos.-16 

Cylinder 

82.9 

76.3 

0.0121 

0.0131 

Viking 

Octagon 

128 

30.8 

0.0078 

0.0325 

Explorer- 11 

Octagon 

203 

72.6 

0.0049 

0.0138 

Explorer- 17 

Sphere 

152 

152 

0.0066 

0.0066 

Sp.  Teles. 

Cylinder 

192 

29.5 

0.0052 

0.0339 

OSO-7 

9-sided 

437 

165 

0.0023 

0.0061 

OSO-8 

Cylinder 

147 

47.2 

0.0068 

0.0212 

Pegasus- 3 

Cylinder 

181 

12.1 

0.0055 

0.0826 

Landsat-1 

Cylinder 

123 

25.2 

0.0081 

0.0397 

ERS-1 

Box 

135 

12.0 

0.0074 

0.0833 

LDEF-1 

12-face 

169 

93.1 

0.0059 

0.0107 

HEAO-2 

Hexagon 

174 

80.1 

0.0057 

0.0125 

Vanguard- 2 

Sphere 

23.5 

23.5 

0.0426 

0.0426 

SkyLab 

Cylinder 

410 

47.1 

0.0024 

0.0212 

Echo-1 

Sphere 

0.515 

0.515 

1.9417 

1.9417 
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Per  Table  2,  B*  =  0.01  will  be  chosen  as  a  nominal  value  for  this  research. 

The  location  of  the  perigee  and  apogee  were  calculated  using  the  state  integrated 
over  6000  time  units  (TU).  Visible  in  Figure  3,  perigee  and  apogee  exhibit  the  expected 
secular  behavior.  While  not  proof  that  air  drag  effects  were  correctly  incorporated 
into  the  equations  of  motion,  the  plots  offer  a  quick  common-sense  check  because  the 
apogee  decreases  and  the  perigee  remains  relatively  constant  over  time. 


+1.152 


=  0.0025 
O 


30 

Time  [days] 


Figure  3:  Apogee  and  Perigee  Location  Over  Time  in  the  Prescence  of  Air  Drag 


3. 2. 2. 3  Third-Body  Moon.  The  motion  of  the  Moon  is  indeed  complex. 
Accurate  modeling  of  the  Moon’s  orbit  is  complicated  by  third-body  gravity  from  the 
Sun,  as  well  as  the  Moon’s  inclination  of  5.1°  from  the  Earth’s  equator.  The  Moon’s 
motion  can  be  simplified  by  assuming  that  its  orbital  plane  is  aligned  with  the  Earth’s 
equator  and  by  assuming  that  the  Moon’s  orbit  is  circular.  These  assumptions  are  valid 
because  the  primary  intent  of  this  research  is  to  gain  understanding  of  how  the  Moon 
alters  the  motion  of  a  satellite  near  the  KAM  Torus.  In  future  applications  with  the 
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intent  of  accurately  predicting  satellite  states,  a  more  complex  Moon  model  will  be 
required. 

Considering  the  Moon’s  position  is  required  in  the  ECEF  frame,  it  can  be  as¬ 
sumed  that  the  Moon  begins  positioned  on  the  ECEF  x-axis  at  60.27  Earth  radii 
(60.27  distance  units  (DU))  and  follows  a  circular  path  about  the  Earth,  as  is  seen  in 
Figure  4. 


Figure  4:  Earth,  Moon,  and  ECI/ECEF  Coordinate  Frames 

Care  must  be  taken  to  correctly  determine  the  Moon’s  orbital  frequency  with 
respect  to  the  ECEF  frame.  Because  the  Earth’s  rotational  frequency  is  larger  than 
the  Moon’s  orbital  frequency,  the  differential  between  the  two,  0Jmoon  —  a;©,  results  in 
the  Moon  rotating  clockwise  about  the  Earth  in  the  ECEF  frame: 


ECEF 

moon 


T'moon 


COS  \{tdrnoon 
sill  \{i^rnoon 


(U©)t] 

ue)t] 


0 


(52) 


where  ujmoon  is  the  Moon’s  position  vector  in  the  earth-centered-inertial  (ECI)  frame, 
n;©  is  the  Earth’s  rotational  angular  velocity  in  the  ECI  frame,  rmoon  is  the  magnitude 
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of  the  Moon’s  position  vector,  and  ECEFrmoon  is  the  Moon’s  position  vector  in  the 
ECEF  frame. 

Next,  the  vector  from  the  satellite  to  the  Moon  must  be  calculated.  Referring 
to  Figure  5,  simple  vector  addition  and  subtraction  yield  an  expression,  Equation  53, 
for  the  vector  between  the  satellite  and  the  Moon. 


Figure  5:  Determining  the  Position  Vector  from  the  Satellite  to  the  Moon 


^  sat /moon  ^ moon  ^ sat 


(53) 


Referring  to  Equation  2  and  applying  F  =  ma,  the  acceleration  of  the  satellite 
due  to  the  Moon  becomes: 


^■s  at  /  moon 


l^moon^ sat /  moon 
rp  3 

7  sat /moon 


(54) 


The  acceleration  can  then  be  added  to  the  Hamiltonian  partials  shown  in  Equa¬ 
tions  48,  49,  and  50. 
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3.3  Torus- Perturbed  Motion  Comparison  Loop 


In  the  comparison  loop,  two  initially  aligned  states  are  updated  and  followed  in 
time.  The  reference  torus  state  is  propagated  by  simply  updating  the  torus  coordinates 
with  the  frequency-time  relationship  described  by  Equation  29.  The  perturbed  state  is 
propagated  via  numerical  integration  as  described  in  §§3.2. 1-2.  Again,  it  is  important 
to  note  that  both  the  KAM  Torus  and  the  numerical  integration  include  the  Earth’s 
geopotential.  The  only  perturbations  considered  in  this  research  are  air  drag  and 
moon  third-body  effects.  During  each  iteration  of  the  comparison  loop,  the  difference 
between  the  states  in  both  physical  and  torus  variables  5q,  5p,  hQ,  and  hP,  are 
calculated  and  recorded  in  a  text  file.  Further  discussion  in  this  research  will  consider 
the  following  two  definitions: 


SX  = 


5q 

5p 


(55) 


hY  = 


<5Q 


SP 


(56) 


3.3.1  Collecting  the  Physical  State  Differences.  Since  both  states  have  been 
updated  after  some  time,  At,  the  difference  between  the  states  in  physical  variables 
can  be  determined.  Currently  the  torus  state  is  known  in  torus  variables  and  the 
perturbed  state  is  known  in  physical  variables.  To  obtain  the  torus  state  in  physical 
variables,  the  same  procedure  utilized  in  §3.2.1  to  extract  the  initial  state  from  the 
torus  is  followed.  Now  with  both  states  described  in  physical  variables,  the  state 
differences  can  be  calculated,  expressed  as  shown  in  Equation  57: 
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q 

p 

integration 

where  q  are  the  physical  coordinates  and  p  are  the  physical  momenta. 

3.3.2  Transforming  Between  Physical  and  Torus  Variables.  To  understand 
the  perturbation’s  effects  on  the  torus  coordinates  and  momenta,  the  data  described 
in  the  physical  variables  must  be  transformed  back  to  the  torus  variables  and  vice 
versa.  This  is  accomplished  with  the  Jacobian  Matrix,  |^.  An  estimate  for  this  matrix 
can  be  obtained  by  defining  three  expressions  of  the  state: 
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(60) 
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where  X  is  the  physical  state  vector,  Y  is  the  canonical  state  vector,  and  Z  is  the 
classical  orbital  elements  state  vector. 


Then  it  follows  that  the  Jacobian,  the  linear  transormation  between  the  X  and 
Y  states  can  be  represented  as  Equation  61: 


dY  _  dY  dZ 
dX  ~  ~dZdX 

Wiesel  illustrates  the  development  of  these  Jacobians  in  his  work  discussing  the 
motion  of  a  satellite  near  a  torus  (19).  The  transformations  are  estimates  and  cause 
some  error  when  comparing  against  reality.  This  error  proceeds  from  the  fact  that  the 
state  vectors  are  described  with  only  the  second  zonal  harmonic  of  the  geopotential.  To 
provide  the  Y  state  vector  with  utmost  accuracy,  the  entire  geopotential  would  have 
to  be  considered  when  developing  the  Jacobian  which  is  impossible  to  do  analytically. 

The  state  differences  in  terms  of  the  torus  coordinates  can  be  estimated  as  in 
Equation  62: 


1 

G* 

1 _ 

dY 

Jq 

SP 

~  dX 

Sp 

The  state  differences,  SQ  and  JP,  are  now  available  for  analysis. 


(62) 
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3.3.2. 1  Improving  the  Accuracy  of  the  Jacobian.  In  an  effort  to  make 
more  accurate,  it  is  observed  that  its  inverse,  can  be  expanded  to  Equation  63: 


<9X 

dY 


dq  dq 
9Q  OP 

dp  dp 
dQ  OP 


where  it  is  seen  that  the  left-hand  side  of  the  matrix 
Fourier  series.  At  this  time  the  author  is  unaware  of 
exact  values  for  the  right-hand  side  of  the  matrix  in 
available,  it  would  certainly  be  utilized. 


(63) 

is  available  from  the  torus  via  the 
a  practical  approach  to  obtaining 
Equation  63.  If  an  approach  were 


The  left  hand  portion  of  the  matrix  described  in  Equation  63  is  presented  for 
both  cases  at  the  initial  time,  t  —  0,  in  Table  3. 


Table  3:  Comparison  of  Portion  of  Jacobian  Approximated  by  Two-Body  Problem 
and  Exact  Values  from  Torus  ffl 


Two-Body  Jacobian 

Partial  Torus  Extraction  Jacobian 

-6.99964e-l 

-7.33130e-l 

-6.67780e-l 

-6.99133e-l 

-7.33187e-l 

-6.67433e-l 

-7.85493e-l 

-6.91869e-l 

-8.24930e-l 

-7.84273e-l 

-6.91947e-l 

-8.23960e-l 

-3.34563e-l 

-1.71706e-15 

-3.12896e-l 

-3.33773e-l 

-4.24398e-5 

-3.12385e-l 

5.93428e-l 

6.75994e-l 

5.64685e-l 

5.92915e-l 

6.76079e-l 

5.64441e-l 

-6.28819e-l 

-6.02387e-l 

-6.64143e-l 

-6.28013e-l 

-6.02444e-l 

-6.63369e-l 

3.91351e-l 

1.00684e-15 

3.77868e-l 

3.91511e-l 

-1.12216e-5 

3.78154e-l 

A  comparable  presentation  of  data  is  available  for  the  ISS  Torus  and  is  located 
in  Appendix  A,  Table  11.  In  both  cases,  significant  updates  are  made  to  the  Jacobian 
in  some  places  and  in  other  places  the  two-body  estimate  is  very  close  to  the  actual 
value  provided  by  the  torus.  Once  the  update  to  the  left-hand  side  of  the  Jacobian 
has  been  obtained,  the  inverse  must  be  taken  again  as  is  the  Jacobian  required  for 
the  linear  transformation  from  physical  variables  to  torus  variables.  Throughout  the 
remainder  of  this  work,  the  two  versions  of  the  Jacobian  will  be  referred  to  as  the 
partial  torus  extraction  (PTE)  Jacobian  and  the  two-body  problem  (2BP)  Jacobian. 
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3. 3. 2. 2  Illustration  of  Jacobian  Accuracy  Effects.  For  the  sake  of  com¬ 
paring  the  two  Jacobian  versions,  a  relationship  between  the  derivative  of  the  physical 
state  and  the  derivative  of  the  torus  state  is  derived.  Starting  with  the  linear  trans¬ 
formation  between  state  differences 


SY  =  Msx  (64) 

where  5Y  is  the  difference  vector  in  torus  variables  between  the  perturbed,  integrated 
state  and  the  torus  state,  and  5X  is  the  difference  vector  in  physical  variables  between 
the  perturbed,  integrated  state  and  the  torus  state. 

Next,  since  the  state  derivative  relationship  is  desired,  the  time  derivative  is 
applied  to  Equation  64: 


JY 


d_ 

dt 


ax 

ax'® 


(65) 


Now  that  the  chain-rule  has  been  applied  it  is  apparent  that  the  state  derivative 
is  dependent  on  the  time  derivative  of  the  Jacobian: 


dY  ■ 

5X+—5X 

oX 


Selectively  expanding  the  differences  yields 


(66) 


int 


V  -  -  (—\  AY 

torus  ~  dt\dx  r  +  ox 


xint  -  x 


torus 


(67) 


where  Yint  is  the  perturbed,  integrated  state  derivative  transformed  into  torus  state 
variables,  Xint  is  the  perturbed,  integrated  state  derivative  in  physical  state  variables, 
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Y torus  is  the  reference  torus  state  derivative  transformed  into  torus  state  variables, 
and  Xi:orus  is  the  reference  torus  state  derivative  in  physical  state  variables. 

Through  selective  application  of  the  definition  of  the  time  derivative: 


Y  --Y  -  - 

A  int  o .  ±  torus  7 , 

at  at 


fiX  + 


dY 

BX 


B 

Bt 


X 


torus 


Next,  rearranging 


(68) 


V  -  d  (dY\w+dYii  dYdv  dv 

in,~  dt  fax!  +ax^”*  dxdt  ,m‘‘  dt  ,m‘‘ 


With  the  cancellation  of  terms,  the  expression  becomes 


(69) 


d  f  BY  \  BY  ■ 

Y‘”-s(ax)ix+axx“  <70> 

Equation  70  has  an  important  ramification.  When  perturbed  motion  is  inte¬ 
grated  and  compared  against  the  reference  torus  state  some  state  difference  5X  will 
necessarily  exist.  The  state  difference  doesn’t  make  Equation  70  unusable,  but  some 
complexity  is  incurred  in  the  requirement  to  calculate  the  time  derivative  of  the  Ja¬ 
cobian.  With  the  additional  choice  of  evaluating  the  Jacobian  with  the  unperturbed, 
integrated  motion  the  additional  simplification,  SX  =  0,  can  be  made.  This  simplifi¬ 
cation  yields  Equation  71: 


BY  ■ 

v  —  _ _ Y 

int  (9X 


(71) 


where  X  is  available  as  Equations  45,  46,  47,  48,  49,  and  50.  Again,  Equation  71 
is  only  valid  when  the  integration  and  torus  states  (coordinates  and  momenta)  are 
aligned  near-perfectly,  otherwise  Equation  70  must  be  used. 
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The  time  derivatives  of  the  toms  coordinates  of  the  unperturbed  motion  near 
the  reference  torus  were  calculated  utilizing  the  two-body  problem  Jacobian  as  well 
as  the  Jacobian  updated  via  the  torus  extraction.  It  is  expected  that  the  coordi¬ 
nate  derivatives  should  closely  approximate  the  torus  basis  frequencies.  The  results 
shown  in  Figure  6  illustrate  that  the  PTE  Jacobian  yields  more  reasonable  results  for 
instantaneous  Qi  values  than  the  two-body  Jacobian.  The  PTE  Jacobian  result  ex¬ 
hibits  small  oscillations  about  the  basis  frequency,  ui,  while  the  2BP  Jacobian  result 
oscillates  to  a  greater  degree. 


Figure  6:  Time  Derivative  of  Unperturbed  Torus  Coordinate  Q i  Utilizing  the  Two- 
Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 


Even  more  importantly,  the  averaged  value  of  Qi  from  the  PTE  Jacobian  is 
closer  to  the  basis  frequency  than  the  averaged  value  of  Q i  from  the  2BP  Jacobian. 
The  results  for  all  three  coordinate  rates  are  summarized  in  Table  4. 
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Table  4:  Averages  of  Coordinate  Time  Derivatives  Calculated  Using  the  Two-Body 
Jacobian  &  the  Partial  Torus  Extraction  Jacobian 


Qx  [rad/TU] 

Q2  [rad/TU] 

Q3  [rad/TU] 

Torus  Frequency 

0.86115953 

-0.05983301 

0.00158568 

2BP  Jacobian 

0.86123874 

-0.05983382 

0.00150735 

PTE  Jacobian 

0.86117945 

-0.05983281 

0.00156017 

In  addition,  for  the  unperturbed  (no  perturbing  forces  other  than  the  geopoten¬ 
tial)  motion  near  the  reference  torus,  it  is  expected  that  P  =  0.  The  plots  in  Figure  7 
demonstrate  that  while  both  the  2BP  and  PTE  Jacobian  results  oscillate  about  zero 
for  Pi,  the  magnitude  of  the  oscillations  produced  using  the  PTE  Jacobian  are  much 
closer  to  zero. 


Figure  7:  Time  Derivative  of  Unperturbed  Torus  Momentum  PI  Utilizing  the  Two- 
Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 


Comparable  results  are  seen  in  plots  for  Q2,  Q 3,  P2,  and  P3  which  are  included 
in  Appendix  A,  Figure  41,  Figure  42,  Figure  43,  and  Figure  44.  The  average  time 
derivatives  for  the  momenta  using  both  Jacobians  are  summarized  in  Table  5. 
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Table  5:  Averages  of  Momenta  Derivatives  Calculated  Using  the  Two-Body  Jacobian 
&  the  Partial  Torus  Extraction  Jacobian 


Pi  [rad/TU] 

P2  [rad/TU] 

P3  [rad/TU] 

Ideal  Torus  Value 

0 

0 

0 

2BP  Jacobian 

9.3  x  10"7 

8.4  x  10"7 

9.4  x  10"7 

PTE  Jacobian 

1.0  x  lO"6 

8.7  x  10"7 

1.0  x  10"6 

While  the  values  for  time  derivatives  of  the  momenta  are  all  practically  zero, 
the  2BP  Jacobian  case  appears  slightly  more  accurate  than  the  PTE  Jacobian  case. 
All  factors  considered,  the  PTE  Jacobian  is  clearly  preferred  for  two  reasons.  First, 
the  coordinate  time  derivatives  more  closely  approximate  the  torus  basis  frequencies. 
Second,  at  any  given  time  the  values  of  the  coordinates  derivatives  and  momenta 
derivatives  remain  closer  to  the  theoretical  torus  values  described  by  the  basis  fre¬ 
quencies.  As  mentioned  earlier,  while  the  Jacobians  and  were  improved  by 
updating  the  left-hand  side  of  important  accuracy  improvements  may  be  gained 
if  a  method  for  obtaining  exact  values  for  the  right-hand  side  is  developed. 

The  effects  of  the  improvement  of  the  Jacobian  are  even  more  dramatically  seen 
when  applied  to  Torus  #2,  the  ISS  Torus,  as  shown  in  Figure  8. 

The  coordinate  rate  Qi  which  is  the  frequency  of  mean  anomaly,  wildly  oscillates 
between  0  [rad/TU]  and  2.5  [rad/TU]  in  the  two-body  problem  Jacobian  estimate 
case.  The  actual  frequency,  0.92  [rad/TU]  is  much  more  closely  approximated  via  the 
PTE  Jacobian.  Note  that  a  zero  value  for  the  rate  of  change  of  mean  anomaly  infers 
that  the  satellite  is  essentially  stopped  in  it’s  orbit  for  a  short  period  of  time  which  is 
not  realistic  and  tells  us  that  the  two-body  estimate  is  not  accurate  in  comparison  to 
the  PTE  Jacobian.  Similar  results  are  seen  in  plots  for  Q2,  Q3 ,  A,  A,  and  A  which 
are  included  in  Appendix  A,  Figures  45,  46,  47,  48,  and  49. 
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Figure  8:  Time  Derivative  of  Unperturbed  Torus  #2  (ISS)  Coordinate  Qi  Utilizing 
the  Two-Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 

3.3.3  Update  of  the  Torus  Phase.  Nominally,  the  coordinates  of  the  KAM 
Torus  update  linearly  with  time,  but  it  must  be  considered  that  motion  near  the  torus 
is  being  evaluated.  When  describing  perturbed  motion  near  the  torus,  the  coordinates 
no  longer  increase  in  an  exclusively  linear  manner. 

Over  time  the  perturbed  motion  will  exhibit  a  non-negligible  phase  difference  in 


each  of  the  coordinates.  The  growth  of  the  magnitude  of  position  error  over  time  for 
unperturbed  motion,  drag  perturbed  motion,  and  third-body  lunar  perturbed  motion 
near  the  reference  torus  is  shown  in  Figure  9. 


Figure  9:  Actual  Error  for  Three  Cases  of  Motion  Near  a  Reference  Torus 


The  error  seen  in  Figure  9  is  not  solely  due  to  phase  offset  but  also  includes  error 
due  to  momenta  differences  between  the  reference  torus  and  the  perturbed  motion. 
The  baseline  case  shows  little  error  accumulation  but  both  perturbed  cases  produce 
substantial  error  over  several  days.  The  error  magnitude  was  calculated  per  Equa¬ 
tion  72  and  Equation  73: 


£q  =  q int  -  q *  orus  (72) 

r error  =  \j T  <^2  (73) 

where  3q  is  the  difference  vector  of  the  physical  coodinates  between  the  perturbed, 
integrated  state  and  the  reference  torus  state,  q int  is  the  physical  coordinate  vector  of 
the  perturbed,  integrated  state,  q torus  is  the  physical  coordinate  vector  of  the  reference 
torus  state,  rerror  is  the  magnitude  of  the  position  error  between  the  two  states,  and  Sx, 
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5y,  and  5z  are  differences  in  physical  coordinates  between  the  perturbed,  integrated 
state  and  the  reference  torus  state. 

Phase  differences  have  an  additional  consequence  beyond  that  described  above. 
As  the  perturbed  state  diverges  from  the  reference  torus  state,  errors  will  be  mag¬ 
nified  because  the  Jacobian  method  used  to  transform  between  physical  and  torus 
coordinates  is  a  linear  transformation.  Over  time,  the  magnitude  of  JQ  at  each  time 
step  will  grow  large  enough  that  the  linear  transformation  using  the  Jacobian, 
will  no  longer  be  valid. 

If  the  phase  discrepancy  is  ignored  and  the  linear  assumption  is  violated  the 
numerical  results  become  invalid.  The  total  position  error  magnitude  over  time  cal¬ 
culated  per  Equations  73  and  74. 


SX 


<9X 

dY 


8  Q 
8P 


(74) 


The  total  position  error  magnitude  over  time  is  shown  in  Figure  10. 


Figure  10:  Error  Magnitude  Calculated  with  Linear  Assumption  Violated 
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The  linear  transformation  barely  holds  for  the  baseline  case  which  means  that 
even  a  very  small  amount  of  discrepancy  between  the  motion  near  the  torus  and  the 
reference  torus  causes  the  linearization  to  magnify  the  error  due  to  differences  in  the 
coordinates.  The  errors  in  the  perturbed  cases  grow  unreasonable  levels  in  just  days 
which  is  not  expected  for  air  drag  or  lunar  third-body  effects  (compare  with  accurate 
error  plots  in  Figure  9).  Per  this  observation,  the  torus  coordinates  absolutely  must  be 
updated  at  each  time  step  if  knowledge  is  to  be  gained  from  analyzing  differences  in 
the  torus  coordinates  and  momenta.  The  torus  coordinates  will  be  updated  as  follows 

Q,torus(t)  QtorusO  T  T  (75) 

where  Q torus  is  the  coordinate  vector  at  a  given  time  t,  Qtoruso  is  the  initial  coordinate 
vector  at  initial  time  to,  is  the  torus  basis  frequency  vector,  and  is  the  phase 
vector  at  a  given  time  t. 

The  phase  vector  is  accumulated  over  time  in  order  to  compare  the  perturbed 
motion  with  the  closest  location  on  the  reference  torus  at  any  time  such  that: 

^  Q  int  Q  torus  (76) 

An  important  distinction  must  be  made  now  that  the  reference  torus  coordinates 
are  updated  to  minimize  error  with  the  perturbed  motion  at  every  time  step.  At  each 
time  step,  additional  phase,  designated  as  5#,  will  be  added  to  the  current  phase  to 
keep  the  two  states,  perturbed  and  reference  torus,  aligned. 

Said  another  way,  is  the  unique  amount  of  phase  added  to  the  total  phase 
at  each  time  step  such  that 
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^  d*  Q  int  Qtorus  ^ 


(77) 


Note  that  while  phase  accumulation,  <h,  helps  to  minimize  the  error  at  any  time 
step,  it  does  not  remove  all  error  as  ciP  has  not  been  accounted  for.  A  difference  in 
torus  momenta,  AP ,  implies  that  the  motion  has  moved  from  the  torus  to  an  adjacent 
torus.  In  this  case,  since  the  motion  near  the  reference  torus  is  under  consideration, 
adjacent  tori  are  not  considered,  but  rather  7P  will  be  evaluated  as  it  pertains  to 
inducing  state  error. 
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IV.  Results 


4-1  Understanding  Differentials  in  the  Torus  Coordinates  and  Momenta 

A  torus  has  coordinates  that  increment  linearly  with  time  and  momenta  that 
remain  constant  with  time.  When  motion  near  a  torus  is  perturbed  by  drag  or  third- 
body  lunar  effects,  non-linear  changes  are  seen  in  the  coordinates  as  additional  phase 
not  present  in  the  reference  torus. 

Per  Equation  77,  at  each  time  step  phase  differences  for  each  coordinate  will  be 
collected  as 

S*Q1 

hd><23 

Clearly  visible  in  Figure  11,  the  Q i  coordinate  grows  faster  over  time  in  the 
drag  case  and  in  the  lunar  case  exhibits  some  type  of  non-linear  behavior. 
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Figure  11:  Coordinate  Q1  Phase  Differences  for  Baseline,  Drag,  and  Lunar  Cases 
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While  Figure  11  makes  it  fairly  obvious  that  the  drag  and  lunar  cases  show  non¬ 
linear  changes  in  the  coordinates,  the  plot  does  not  make  clear  how  the  coordinate 
changes  in  the  long  run.  This  is  remedied  by  observing  the  phase  accumulation  -  the 
amount  of  phase  required  to  most  closely  align  the  state  of  the  reference  torus  with 
that  of  the  state  of  the  perturbed  motion.  The  phase  accumulation  for  Q i  is  shown 
in  Figure  12. 


Figure  12:  Coordinate  Q1  Phase  Accumulation  for  Baseline,  Drag,  and  Lunar  Cases 


In  theory,  the  phase  accumulation  for  the  baseline  case  should  be  zero  for  all 
time,  but  the  plot  has  a  small  level  of  drift  error.  The  source  of  this  error  is  most 
likely  that  the  torus  was  not  an  exact  fit  to  the  integrated  unperturbed  data.  For  the 
drag  case,  the  phase  accumulation  has  a  quadratic  nature  which  is  likely  associated 
with  the  exponential  density  characteristics  of  the  Earth’s  atmosphere.  Since  Q i  is 
approximately  the  mean  anomaly  of  an  orbit,  it  is  reasonable  that  as  an  orbiting 
satellite  loses  energy  due  to  drag,  its  semi-major  axis  is  decreased  and  its  velocity  will 
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increase,  thus  increasing  M.  Again,  the  lunar  case  exhibits  a  deviation  from  linear 
behavior  but  it  is  difficult  to  characterize  at  this  time  scale. 

Referring  back  to  §2. 1.2.1  on  drag  effects,  Equation  20,  it  is  predicted  that  drag 
will  not  affect  the  RAAN.  Q2  is  approximately  equivalent  to  RAAN,  and  in  Figure  13 
it  is  seen  that  additional  phase  is  accumulated  due  to  drag  but  the  magnitude  of  $<52 
is  very  small  compared  to  Q2. 


2  3 

Time  [days] 


Figure  13:  Coordinate  Q2  Phase  Accumulation  for  Baseline,  Drag,  and  Lunar  Cases 


Q 3  which  is  synonymous  with  AP  exhibits  quadratic  behavior  and  again  the 
lunar  case  has  definite  effects  as  well,  as  shown  in  Figure  14. 

Thus  far,  torus  coordinates  have  been  discussed,  but  the  momenta  must  be 
considered  as  well.  The  non-perturbed  torus  should  have  constant  momenta.  Non¬ 
constant  momenta  imply  that  the  satellite  has  moved  from  the  reference  torus  to  an 
adjacent  torus.  An  adjacent  torus  is  characterized  as  similar  to  the  reference  torus 
but  with  slightly  different  frequencies  -  which  means  that  the  torus  coordinates  will 
increment  linearly  at  different  basis  frequencies. 
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Figure  14:  Coordinate  Q3  Phase  Accumulation  for  Baseline,  Drag,  and  Lunar  Cases 


Figure  15  shows  that  as  expected,  the  Pi  momentum  is  fairly  constant  for  the 
baseline  case,  while  the  drag  and  lunar  cases  exhibit  a  linear  rate  of  change.  The 
change  in  momenta,  SP2  and  dP3,  behave  similarly  and  are  shown  as  Figure  52  and 
Figure  53  in  Appendix  A. 

Assuming  the  reference  torus  is  used  to  predict  the  location  of  a  perturbed 
satellite  at  a  given  time,  some  state  error  exists  due  to  the  coordinate  and  momenta 
offsets  from  the  reference  torus,  /)<&  and  dP  respectively.  As  detailed  in  §3.3.1  /  Equa¬ 
tion  57,  dq  and  dp  are  available,  so  the  magnitude  of  the  error  can  be  calculated  as 
per  Equation  79  and  Equation  80: 


Terror  =  \/  dx2  +  d?/2  +  dz2 


(79) 


ve„„  =  +  <siv  +  W? 


(80) 


where  rerror  and  Verror  are  the  position  and  velocity  error  magnitudes,  respectively. 
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Figure  15:  Momentum  Differences,  5  Pi,  for  Baseline,  Drag,  and  Lunar  Cases 
The  total  position  and  velocity  error  for  the  three  cases  are  shown  in  Figure  16 


and  Figure  17,  respectively. 


Figure  16:  Total  Position  Error  for  Baseline,  Drag,  and  Lunar  Cases 
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Figure  17:  Total  Velocity  Error  for  Baseline,  Drag,  and  Lunar  Cases 


The  oscillations  seen  in  the  drag  case  in  Figure  16  approach  zero  kilometers 
error  in  every  cycle.  This  phenomenon  is  associated  with  the  satellite  being  located 
at  the  perigee.  The  errors  for  the  drag  and  lunar  cases  seem  large,  but  it  is  important 
to  note  that  3-10  km  of  error  when  the  position  vector  magnitude  is  on  the  order 
of  7000  km  is  just  over  1%  error.  It  is  tempting  to  believe  that  this  error  can  be 
reduced  by  accurately  predicting  the  effects  of  drag  and  lunar  third-body  on  the 
torus.  The  approach  followed  thus  far  has  already  updated  the  torus  coordinates  after 
each  time  step.  Therefore,  the  remaining  error  is  present  in  the  momenta  and  cannot 
be  removed  with  a  state  update  without  jumping  to  an  adjacent  torus.  Rather,  this 
approach  has  demonstrated  that  updating  the  torus  coordinates  is  not  just  beneficial 
but  absolutely  necessary  in  applications  that  require  use  of  the  Jacobian,  and  that 
the  remaining  error  resides  in  the  momenta  differences. 


48 


4-2  Isolating  Primarily  Effected  Torus  Variables 

To  understand  which,  if  any,  of  the  specific  coordinates/momenta  are  most  af¬ 
fected  by  the  perturbations,  assume  a  change  in  only  one  coordinate  or  momentum 
and  calculate  the  associated  error.  For  example: 


™-  ax 

SX=SY 


0 

0 

0 

0 

0 


(81) 


Note  that  it  has  already  been  shown  that  differences  in  the  coordinates  cause 
large  errors  and  therefore  the  torus  coordinates  are  updated  with  phase;  therefore,  it 
is  expected  that  the  coordinates  should  not  contribute  much  error  since  at  any  given 
time  step  the  only  error  in  the  coordinates  is  <5$  as  $  is  already  accounted  for.  The 
total  error  per  associated  coordinate  or  momentum  is  calculated  using  Equations  79 
and  80  .  Figures  18,  19,  20,  21,  22,  and  23  show  the  total  position  error  over  time 
due  to  SPi,  5P2,  and  5P3,  respectively. 
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Figure  18:  Total  Position  Error  Due  to  for  Baseline,  Drag,  and  Lunar  Cases 


Figure  19:  Total  Position  Error  Due  to  dd>Q2  for  Baseline,  Drag,  and  Lunar  Cases 
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Figure  20:  Total  Position  Error  Due  to  ^$<53  for  Baseline,  Drag,  and  Lunar  Cases 


Figure  21:  Total  Position  Error  Due  to  8P\  for  Baseline,  Drag,  and  Lunar  Cases 
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Figure  23:  Total  Position  Error  Due  to  SP3  for  Baseline,  Drag,  and  Lunar  Cases 
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As  expected,  the  error  resulting  from  each  coordinate  is  minimal,  mainly  due 
to  the  approach  of  updating  the  coordinate  phase  after  each  time  step.  The  error 
contributed  from  SP\  and  SPo  is  on  the  same  order  as  the  coordinates,  but  the  error 
due  to  SP3  is  significant.  This  indirectly  confirms  the  observation  that  updating  the 
coordinates  is  necessary  and  error  resides  in  the  momenta  offsets. 

4-3  Methods  to  Characterize  Air  Drag  Effects  on  KAM  Torus  Coordi¬ 
nates  and  Momenta 

It  is  not  sufficient  to  be  able  to  numerically  integrate  a  specific  case  or  even  to 
view  data  collected  in  the  past  to  see  how  a  satellite  has  traversed  or  deviated  from 
a  reference  torus.  The  capability  sought  is  to  be  able  to  predict  the  effects  of  air  drag 
and  third-body  lunar  effects  on  motion  near  a  given  torus. 

With  that  thought  in  mind,  it  is  intuitive  to  attempt  to  establish  a  method 
of  coordinate  rate  estimating  by  manipulating  the  results  discussed  in  §4.1.  Two 
approaches,  based  on  the  same  data  set,  will  be  considered.  First,  the  coordinate 
differences  may  be  divided  by  the  time-step  and  then  fit  via  linear  least  squares 
resulting  in  coordinate  time  derivatives  that  are  functions  of  time.  Second,  the  phase 
accumulation  for  each  coordinate  can  be  divided  by  the  associated  time  vector  then 
fit  via  linear  least  squares  to  provide  a  function  returning  the  phase  estimate  for  each 
coordinate  at  a  given  time. 

4-3.1  Estimated  Differences  of  Coordinate  Time  Derivatives.  Time  deriva¬ 
tives  of  the  coordinates  are  useful  for  applications  that  require  numerical  integration 
and  also  for  comparison  against  the  basis  frequencies  of  the  reference  torus.  In  the 
comparison  loop  described  in  §3.3,  the  phases  of  the  coordinates  are  updated  at  every 
time  step.  This  update  results  in  a  collection  not  of  the  difference  between  unper- 


53 


turbed  and  perturbed  coordinates,  but  rather  a  phase  difference,  5#,  between  the 
current  and  previous  time  step  that  does  not  include  the  summed  phase  from  all  pre¬ 
vious  time  steps.  The  differences  between  the  reference  torus  coordinate  derivatives 
and  the  perturbed  motion  coordinate  derivatives  are  calculated  as 

**  =  M  <82> 

The  results  of  approximating  the  coordinate  time  derivative  differences  for  two 
different  time  periods,  500  TU  and  6000  TU  (about  5.5  days  and  56  days),  are  sum¬ 
marized  in  Table  6. 


Table  6:  Intercept  and  Slope  for  Torus  #1  Coordinates  Rates  Via  Linear  Least  Squares 


500  TU@  At  =  0.25  TU 

6000  TU  @  At  =  0.25  TU 

Coordinate 

Intercept 

Slope 

Intercept 

Slope 

Rate 

[rad/TU] 

[rad/TU2] 

[rad/TU] 

[rad/TU2] 

dQi 

-3.006  x  10“6 

4.782  x  10“7 

-2.440  x  10“5 

4.800  x  10“7 

sq2 

1.322  x  10“8 

-1.250  x  10"9 

6.831  x  10“8 

-1.244  x  10"9 

6Q3 

5.153  x  10“7 

-1.884  x  10“10 

-1.182  x  10“7 

1.995  x  10“9 

Comparing  the  results  for  the  linear  least  squares  fits  for  the  time  periods  of 
500  TU  and  6000  TU,  it  is  seen  that  the  different  time  spans  don’t  make  a  signif¬ 
icant  difference  especially  with  respect  to  the  function  slopes  of  SQ\  and  SQo-  The 
intercept  values  seem  to  vary  but  they  are  on  the  order  of  the  slopes  and  are  there¬ 
fore  insignificant  after  even  a  few  time-steps.  The  time  derivative  of  coordinate  SQ3 
exhibits  behavior  different  from  the  other  two  -  unexpectedly,  the  slope  for  the  500 
TU  integration  is  negative  while  the  slope  for  the  6000  TU  integration  is  positive. 
The  raw  data  with  the  linear  fits  superimposed  for  all  three  coordinates  for  both  500 
TU/4.5  days  and  6000  TU/56  days  are  shown  in  Figure  24  and  Figure  25. 

Of  particular  interest  is  the  long  period  oscillation  in  Q3  beginning  to  be  visible 
in  the  56  day  plot.  The  non-symmetrical  nature  of  oscillation  certainly  explains  why 
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Figure  24:  Coordinate  Time  Derivative  Differences  and  Linear  Fits  for  4.5  Days 
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Figure  25:  Coordinate  Time  Derivative  Differences  and  Linear  Fits  for  56  Days 
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the  linear  slope  of  the  least-squares  fit  is  changing  from  negative  to  positive  in  the  4.5 
day  and  56  day  cases,  respectively. 

With  the  linear  least  squares  fit  results,  the  torus  coordinate  rates  can  now  be 
described  as 


Q( 0  &  ^Qinterceptt  4“  sloped  (83) 

where  5Qintercept  is  the  linear  least  squares  intercept  for  5Q  and  SQsiope  is  the  linear 
least  squares  slope  for  <5Q. 

Equation  83  can  now  be  used  for  applications  in  which  numerical  integration  of 
the  torus  coordinates  is  required. 


4-3.2  Coordinate  Phases  at  a  Given  Time.  During  the  comparison  loop  de¬ 
scribed  in  §3.3  the  phase  history  previously  termed  phase  accumulation,  was  collected 
over  time.  This  phase  history  appears  to  have  a  quadratic  nature  when  plotted.  When 
the  phase  time  history  is  divided  by  the  time  vector,  the  resulting  data  is  linear  and 
appropriate  for  linear  least  squares.  The  slope  and  intercept  values  from  the  linear 
least  squares  fit  are  summarized  in  Table  7  and  the  approximate  functions  are  com¬ 
pared  against  the  raw  data  in  Figure  26  and  Figure  27. 

Table  7:  Intercept  and  Slope  for  Torus  #1  Coordinates  Phase  Per  Time  Via  Linear 
Least  Squares 


500  TU@  At  =  0.25  TU 

6000  TU  @  At  =  0.25  TU 

Coordinate 

Phase 

Intercept 

[rad/TU] 

Slope 

[rad/TU2] 

Intercept 

[rad/TU] 

Slope 

[rad/TU2] 

&Ql/t 

-4.626  x  10“6 

2.441  x  10“7 

-7.646  x  10“6 

2.358  x  10“7 

^Ql/t 

-6.977  x  10“9 

-5.649  x  10“10 

1.976  x  10“8 

-6.098  x  10”10 

^Qs/t 

1.644  x  10“6 

-3.478  x  10“9 

2.142  x  10“7 

9.147  x  10“10 
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Figure  26:  Coordinate  Phase  Per  Time  and  Linear  Fits  for  4.5  Days 
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Figure  27:  Coordinate  Phase  Per  Time  and  Linear  Fits  for  56  Days 
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The  intercept  and  slope  information  can  be  used  to  incorporate  air  drag  ef¬ 
fects  directly  into  KAM  Torus  applications  for  this  specific  satellite,  orbit,  and  time, 
without  the  need  for  numerical  integration  per  Equation  84: 

Q(t)  A  intercept  A  4? sloped)  t  (84) 

where  $ intercept  is  the  linear  least  squares  intercept  for  and  $>sioPe  is  the  linear  least 
squares  slope  for  <f>. 

Various  trade  studies  and  analyses  would  have  to  be  performed  to  determine 
how  long  the  combined  KAM  Torus  -  Air  Drag  Characterization  would  provide  valid 
solutions.  As  discussed  earlier,  this  coordinate  update  does  not  account  for  differences 
in  momenta  which  do  contribute  considerably  to  error  in  state  knowledge. 

4-3.3  Momenta  Differences.  The  momenta  data  collected  from  the  compar¬ 
ison  loop  differs  from  the  coordinate  data  in  that  it  is  not  utilized  at  any  time  to 
update  the  reference  torus  state.  Because  the  momenta  difference  data  is  approxi¬ 
mately  linear,  it  can  be  fit  via  linear  least  squares  to  provide  a  method  to  predict 
the  momenta  difference  between  the  perturbed  motion  and  the  reference  torus  at  any 
time.  The  function  output  by  linear  least  squares  is  summarized  in  Table  8. 

Table  8:  Intercept  and  Slope  for  Torus  ffl  Momenta  Differences  at  a  Given  Time  Via 
Linear  Least  Squares 


500  TU@  At  =  0.25  TU 

6000  TU  @  At  =  0.25  TU 

Momentum 

Difference 

Intercept 

Slope 

Intercept 

Slope 

<5Pi 

9.805  x  10"7 

-1.935  x  10-7 

9.778  x  10"6 

-1.950  x  10"7 

<5P2 

8.300  x  10“7 

-1.507  x  10“7 

8.074  x  10“6 

-1.527  x  10“7 

SP3 

9.511  x  10“7 

-1.766  x  10“7 

9.179  x  10“6 

-1.781  x  10“7 
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In  the  case  of  the  momenta,  the  intercept  and  slope  values  are  fairly  independent 
of  the  integration  duration.  The  momenta  differences  between  the  reference  torus  and 
the  perturbed  motion  can  easily  be  calculated  as 


*5P(^)  ^P intercept  T  hP sloped  (85) 

where  SPintercept  is  the  linear  least  squares  intercept  for  r)P  and  5Psiope  is  the  linear 
least  squares  slope  for  <5P. 

For  all  three  cases,  over  time,  the  momenta  difference  drifts  negatively  and 
linearly.  The  drift  is  clearly  visible  in  Figure  28  and  Figure  29. 


Figure  28:  Momenta  Differences  and  Linear  Fits  for  4.5  Days 


For  long  term  position  error,  this  infers  that  the  state  error  due  to  instantaneous 
momenta  offsets  will  be  unbounded  over  time  if  just  coordinate  updates  are  performed. 
The  decreasing  slopes  are  representative  of  the  fact  that  momenta  are  energy-like 
terms  which  would  be  expected  to  decrease  over  time  in  the  presence  of  the  non¬ 
conservative  force  due  to  air  drag.  These  insights  are  associated  with  the  necessary 
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Figure  29:  Momenta  Differences  and  Linear  Fits  for  56  Days 


conclusion  that  despite  coordinate  updates/prediction,  a  KAM  Torus  can  only  provide 


accurate  state  prediction  for  a  satellite  encountering  drag  for  a  finite  time  -  unless 


momenta  differences,  hP,  even  though  small,  can  be  accounted  for. 


4-3-4  Differences  Due  to  Time  Period  of  Concern.  As  discussed  in  §§4.3.3- 
4,  the  slopes  of  the  linear  least  square  lines  vary  slightly  depending  on  the  time 
duration  of  the  data  that  was  fit.  Additionally,  in  the  case  of  the  coordinate  Q 3,  the 
slope  even  changes  directions  between  the  4.5  day  and  56  day  cases  of  the  coordinate 
time  derivative  differences  and  coordinate  phase  per  time  predictions.  These  slope 
differences  are  caused  mostly  by  long-periodic  behavior  in  the  coordinates.  This  long- 
periodic  behavior  is  clearly  seen  in  Figure  30.  Phase  and  momenta  difference  plots  for 
the  280  day  case  are  located  in  Appendix  A,  Figures  54  and  55,  respectively. 

Short-term  fits  may  place  emphasis  on  a  peak  or  trough  of  the  periodic  behavior. 
As  a  consequence,  when  utilizing  coordinate  update  techniques  it  must  be  determined 
what  period  of  time  is  of  interest.  If  the  KAM  torus  is  to  be  reconstructed  every 
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Figure  30:  Coordinate  Time  Derivative  Differences  and  Linear  Fits  for  280  days 


couple  of  days,  the  near-term  fit  will  provide  better  accuracy.  Otherwise,  if  the  KAM 
torus  is  to  be  constructed  and  utilized  for  long  time  periods  the  long-term  fit  should 
be  used.  If  both  near-term  and  long-term  accuracy  are  required,  it  may  be  possible 
to  schedule  prediction  curves.  The  coordinate  prediction  would  initialize  with  the 
near-term  prediction  line  and  switch  to  the  long-term  prediction  line  at  a  specified, 
predetermined,  elapsed  time. 


4-3.5  Application  to  Any  Earth  Satellite  Torus.  Once  the  coordinates  and 
momenta  have  been  characterized  as  a  function  of  time,  they  can  simply  be  added 
to  the  propagation  loop  in  future  KAM  torus  algorithms.  The  functions  are  unique 
because  they  are  dependent  on  satellite  configuration,  altitude  and  air  density,  and 
orbit  eccentricity.  It  is  possible  that  in  the  future  a  database  could  be  established  that 
provides  coordinate  and  momenta  functions  for  common  orbit  cases. 

The  author  envisions  that  functions  would  be  available  for  orbits  classified  per 
semi-major  axis,  eccentricity,  and  inclination.  The  functions  could  be  normalized  via 
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some  nominal  ballistic  coefficient  then  modified  per  the  actual  ballistic  coefficient.  A 
brief  survey  was  performed  to  determine  the  feasibility  of  such  an  approach. 


The  differences  in  coordinate  rates,  phase  histories,  and  differences  in  momenta 
were  again  fit  via  linear  least  squares  as  in  the  previous  sections,  but  three  additional 
cases  were  evaluated  in  which  the  ballistic  coefficient  value  was  varied.  The  results  of 
the  drag  coefficient  survey  are  summarized  in  Tabie  9. 

Tabie  9:  Torus  Variable  Fit  Slopes  for  Various  Ballistic  Coefficients 


B*  =  0.01 

Fit  Slope 

B*  =  0.02 

Fit  Slope 

B*  =  0.05 

Fit  Slope 

B*  =  0.1 

Fit  Slope 

SQi 

1.195  x  10“7 

2.392  x  10“7 

6.008  x  10“7 

1.213  x  10“6 

<5Q2 

-3.125  x  lO"10 

-6.288  x  lO"10 

-1.587  x  10"9 

-3.190  x  10"9 

sq3 

-4.711  x  10"11 

2.144  x  10"10 

1.288  x  10"9 

2.783  x  10"9 

&Ql/t 

2.441  x  10-7 

4.846  x  10-7 

1.210  x  10"6 

2.431  x  10"6 

&Q2/t 

-5.649  x  lO"10 

-1.189  x  10"9 

-3.072  x  10"9 

-6.247  x  10"9 

&Qs/t 

-3.478  x  10“9 

-2.624  x  10“9 

-2.773  x  10“n 

4.423  x  10“9 

SPi 

-1.935  x  10“7 

-3.877  x  10“7 

-9.747  x  10“7 

-1.968  x  10“6 

5P2 

-1.507  x  10“7 

-3.019  x  10“7 

-7.592  x  10“7 

-1.533  x  10“6 

SP3 

-1.766  x  10“7 

-3.539  x  10“7 

-8.898  x  10“7 

-1.797  x  10“6 

In  summary  of  the  data  in  Table  9,  the  slopes  of  the  linear  fits  for  nearly  all 
parameters  are  related  such  that  the  fit  slope  is  nearly  proportional  to  the  ballistic 
coefficient.  Again,  the  coordinate  Q3  seems  to  be  behaving  in  a  different  manner.  At 
some  point  in  the  ballistic  coefficient  sweep,  the  parameters  associated  with  Q5,  5Q3 
and  <f>Q3/f,  transition  from  a  decreasing  slopes  to  increasing  slopes.  The  cause  is  again 
related  to  short  and  long  term  periodic  behavior. 

The  results  listed  in  Table  9  are  summarized  graphically  in  Figure  31,  Figure  32, 
and  Figure  33  with  the  exception  that  the  linear  fits  of  the  B*  =  0.02  case  were 
removed  for  clarity: 
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Figure  31:  Linear  Fits  of  Coordinate  Time  Derivative  Differences  for  Various  B*  Val¬ 
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Figure  32:  Linear  Fits  of  Coordinate  Phase  for  Various  B*  Values 
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Figure  33:  Linear  Fits  of  Momenta  Differences  for  Various  B*  Values 
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4-4  Methods  to  Characterize  Third- Body  Lunar  Effects  on  KAM  Torus 
Coordinates  and  Momenta 

Previous  data,  plots,  and  discussions  on  third-body  lunar  effects  on  torus  coor¬ 
dinates  and  momenta  have  primarily  focused  on  data  over  a  time  of  500  TU  or  just 
over  4.5  days.  Because  the  Moon’s  orbital  period  about  the  Earth  is  about  27  days  and 
the  force  of  the  Moon’s  gravity  on  the  satellite  is  clearly  a  function  of  the  geometry 
between  the  satellite  and  the  Moon,  it  is  reasonable  to  suspect  that  data  spanning 
one  or  more  moon  period  is  required  to  uncover  the  entirety  of  the  third-body  behav¬ 
ior.  Thus  a  time  period  of  6000  TU  or  about  56  days  was  chosen  for  additional  data 
collection. 

The  changes  in  phase  required  to  keep  the  reference  torus  aligned  with  the  lunar 
perturbed  motion  appears  to  be  symmetrical  but  there  is  non-negligible  accumulation 
of  phase  over  the  long  term.  The  phase  over  time  exhibits  a  periodic-like  behavior  and 
is  illustrated  in  Figure  34: 


Figure  34:  Coordinate  Q1  Phase  Differences  and  Phase  History  for  Lunar  Case 
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It  is  observed  that  after  one  period,  the  phase  doesn’t  return  to  the  starting 
value.  This  discrepancy  is  at  least  partially  caused  by  the  error/offset  shown  in  the 
baseline  case,  shown  in  Figure  35: 


Figure  35:  Comparison  of  Q1  Phase  History  of  Baseline  and  Lunar  Cases 

From  Figure  35  it  can  be  determined  about  half  of  the  offset  in  the  lunar  case 
can  be  explained  by  the  baseline  offset.  Coordinate  Q?,  phase  shows  similar  behavior 
to  that  of  Qi  phase  except  that  the  Q3  phase  returns  to  the  same  value  after  one 
period  as  shown  in  Figure  36. 

The  phase  of  coordinate  Q2  behaves  much  differently  than  the  other  coordinates 
in  that  while  still  showing  periodic  behavior,  the  phase  is  nearly  constantly  decreasing. 
The  phase  difference  and  history  of  Q2  is  shown  in  Figure  37. 
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Figure  36:  Coordinate  Q3  Phase  History  for  Baseline  and  Lunar  Cases 


0.000010 

0.000005 

0.000000 

-0.000005 

-0.000010 

-0.000015 


Figure  37:  Coordinate  Q2  Phase  Differences  and  Phase  History  for  Lunar  Case 
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In  the  case  that  lunar  perturbed  motion  was  being  compared  to  a  reference 
torus  and  the  coordinate  phases  were  not  being  tracked  and  updated  with  time  the 
error  in  Q2  would  grow  unbounded  while  the  error  in  the  other  two  coordinates  would 
significantly  oscillate  and  grow  slowly.  The  momenta  differences  in  the  lunar  case 
show  significant  behavior.  All  three  exhibit  periodic  behavior  with  a  period  of  about 
37  days.  Within  this  time  period,  the  momenta  differences  return  to  zero  twice,  as  seen 
in  Figure  38.  Both  zero  crossings  are  explained  by  symmetry.  When  the  simulation 
of  the  Moon’s  third-body  effects  begins,  some  angular  offset,  less  than  ir  radians, 
between  the  satellite’s  RAAN,  AP,  and  the  Moon’s  position  about  the  Earth  exists. 
The  first  zero  crossing  is  explained  by  this  initial  angle  being  reduced  to  zero  and 
then  duplicated  on  the  “other  side”.  The  second  zero  crossing  is  then  explained  by 
the  Moon-RAAN-AP  system  completing  an  entire  period. 


0.00004 


Figure  38:  Momenta  Differences  for  Lunar  Case 


The  37  day  period  is  also  clearly  seen  in  the  total  position  and  velocity  error 
plots  shown  in  Figure  39.  The  occurrence  of  zero  error  points  aligning  with  the  zero 
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momenta  difference  points  in  Figure  38  is  not  a  surprise.  Because  this  analysis  requires 
the  coordinates  be  updated  to  account  for  phase  differences  (so  that  ultimately  the 
linear  transformation  via  the  Jacobian  is  valid),  the  only  contributor  to  error  is  nec¬ 
essarily  momenta  differences.  When  the  momenta  differences  between  the  perturbed 
motion  and  the  reference  torus  return  to  zero  the  error  also  returns  to  zero.  Obtaining 
this  accuracy,  even  at  two  specific  points  in  time,  is  difficult  in  practice. 


Figure  39:  Total  Position  and  Velocity  Error  Between  Reference  Torus  and  Perturbed 
Motion  for  Lunar  Cases 


Per  Figure  38  and  Figure  39  in  which  the  system  period  appears  to  be  about  37 
days,  the  approximate  frequency  can  be  estimated  as 


2n  2ti 

Uls-est  =  ^p=  37 daj/ .  107.09 TUjday 


=  0.001586 


,rad 

TU 


(86) 


where  u Jis-est  is  the  estimated  lunar  system  frequency,  and  V  is  the  period  of  the 
combined  frequency  system. 
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As  illustrated  in  Figure  40,  the  system  frequency  is  a  function  of  Q2,  Q3,  and 

^ moon  • 


Figure  40:  Earth  Satellite  and  Third-Body  Lunar  Frequencies  in  the  ECEF  Frame 

Per  the  relationship  defined  in  Figure  40,  Equation  87  specifies  the  system  fre¬ 
quencies  observed  in  Figure  38  and  Figure  39  as 


^ Is  Q 2  ^ moon  T  Q3  (87) 

where  cuis  is  the  frequency  of  the  lunar  system. 

The  lunar  system  frequency  is  a  function  of  the  frequency  of  the  second  and 
third  reference  torus  coordinates  as  well  as  the  frequency  of  the  moon,  all  in  the 
ECEF  frame.  The  frequencies  and  periods  are  summarized  in  Table  10: 

Table  10:  Earth  Satellite  Perturbed  by  Lunar  Effects  System  Frequencies 


Q2 

Q3 

^ moon 

Uls 

Frequency 

rad/TU 

-0.059833 

0.001586 

-0.056683 

-0.001564 

Period 

TU 

105.01 

3961.66 

110.85 

4017.38 

The  result  of  4017.38  TU  or  37.5  days  for  the  combined  period  is  equivalent  to 
the  time  picked  directly  off  the  momenta  differences  and  total  error  plots. 
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It  has  been  shown  that  predicting  the  coordinate  phases  is  not  trivial  because 
the  trends  show  periodic  characteristics  but  do  not  return  to  previous  values  in  the 
case  of  Qi  and  Qo-  Also,  it  must  be  remembered  that  the  Moon  model  utilized  in  this 
research  was  simplified.  To  utilize  this  approach  in  operational  applications,  a  more 
accurate  Moon  model  would  need  to  be  used,  and  in  particular  the  initial  state  of  the 
Moon  would  have  to  align  with  reality  rather  than  arbitrarily  established. 

With  the  complex  characteristics  uncovered  regarding  the  Moon’s  effect  on  mo¬ 
tion  near  a  reference  torus,  an  important  realization  is  made:  if  possible  the  application 
of  KAM  Theory  to  Earth  satellites  should  be  expanded  to  include  the  Moon  in  the 
dynamics  when  constructing  a  KAM  torus.  There  are  two  major  reasons  this  approach 
may  be  viable:  1)  the  motion  near  the  reference  torus  is  seen  to  be  affected  by  specific 
frequencies  that  are  functions  of  current  KAM  torus  frequencies  as  well  as  the  Moon’s 
orbital  frequency  (Bordner  (2)  also  saw  frequencies  due  to  third-body  Moon  effects 
in  his  analysis  of  GPS  satellite  tori),  and  2)  the  forces  applied  to  an  Earth  orbiting 
satellite  by  the  moon  are  conservative. 

Beyond  evidence  showing  it  may  be  possible,  there  is  also  significant  motivation 
to  include  the  moon  in  KAM  torus  construction.  First,  this  research  has  shown  that 
especially  in  the  short  term,  lunar  effects  can  cause  more  state  error  than  air  drag  - 
specifically  seen  in  the  case  of  Torus  #1.  The  maximum  state  error  due  to  the  Moon 
is  on  the  order  of  80  km  and  40  m/s  for  the  analyzed  case  which  has  benefited  from 
phase  updates  to  the  coordinates  after  every  time  step. 

The  second  motivation  is  that  to  predict  and  update  the  coordinates  would  be 
a  significant  effort,  unique  to  any  orbit  at  any  time.  In  practice,  the  requirements 
would  essentially  be  on  par  with  special  perturbation  theory,  which  opposes  the  goals 
of  applying  KAM  tori  to  Earth  satellites.  While  complex  to  initially  formulate,  the 
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new  KAM  torus  including  the  Moon’s  dynamics  will  most  likely  be  adjacent  to  the 
previous  KAM  torus  and  will  have  additional  basis  frequency/momentum  pairs. 
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V.  Conclusion 


5. 1  Results 

This  research  has  realized  several  key  results  that  will  help  to  further  the  appli¬ 
cation  of  KAM  Theory  to  Earth  orbiting  satellites.  The  motion  of  a  satellite  perturbed 
by  atmospheric  air  drag  or  third-body  lunar  gravity  effects  was  compared  to  the  mo¬ 
tion  predicted  by  a  reference  orbit  characterized  by  a  KAM  Torus.  It  was  shown  that 
KAM  Theory  isn’t  precluded  in  these  perturbing  situations,  but  certain  considerations 
must  be  made.  Primary  results  include: 

•  Obtaining  an  accurate  Jacobian  is  critical.  The  quality  of  the  Jacobian  can  be 
greatly  enhanced  beyond  that  of  the  2BP  approximation  by  updating  the  first 
three  columns  of  with  data  extracted  directly  from  the  KAM  torus,  and 

.  Methods  should  be  developed  to  more  accurately  approximate  the  last  three 
columns  of  the  Jacobian,  and 

•  The  linear  transformation  performed  via  the  Jacobian  is  only  valid  when  the 
reference  torus  coordinates  are  closely  aligned  with  the  equivalent  coordinates 
of  the  perturbed  motion.  The  coordinate  alignment  can  be  characterized  as  a 
set  of  phases  that  are  updated  over  time. 

•  Given  an  Earth  satellite’s  state  and  drag  characteristics,  functions  of  time  can 
be  numerically  derived  to  approximate  and  predict  air  drag’s  effects  on  the 
reference  torus  coordinates  and  momenta.  Described  in  the  torus  variables,  air 
drag  primarily  affects  the  first  and  third  coordinates  which  are  approximately 
equivalent  to  the  mean  anomaly  and  argument  of  perigee,  respectively.  Momenta 
differences  cause  unbounded  error  over  time  inferring  that  at  some  point  in  time, 
dependent  on  application,  the  current  torus  characterization  will  become  invalid. 

•  Moon  third-body  effects  cause  significant  near-term  state  error.  In  the  long  term, 
some  of  the  Moon’s  effects  average  out  in  a  specific  period  of  time  determined 
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by  the  Moon’s  orbital  frequency  and  the  KAM  Torus’  second  and  third  basis 
frequencies.  Characterizing  the  Moon’s  effects  on  motion  near  a  reference  KAM 
torus  at  any  given  time  will  require  complex  effort  of  similar  specificity  to  that 
of  special  perturbation  theory  if  the  current  method  is  employed. 

5.2  Future  Efforts  and  Recommendations 

Though  several  key  results  have  been  realized,  as  in  most  pursuits  of  advance¬ 
ment,  more  work  needing  attention  has  been  uncovered.  The  motivation  for  this  re¬ 
search  was  to  contribute  to  a  new  method  of  orbit  characterization  and  prediction. 
Current  day  techniques  are  not  adequate  to  predict,  let  alone  prevent,  the  catastrophic 
collisions  of  artificial  satellites  that  will  most  likely  hamper,  if  not  deny,  humanity’s 
access  to  space  in  the  future.  With  consideration  to  the  lofty  goal  of  improving  upon 
the  methods  used  for  decades  to  track  and  predict  Earth  satellite  motion  but  also 
with  the  focus  pertinent  to  this  body  of  research,  the  following  future  efforts  are 
recommended. 

•  Develop  a  database  of  coordinate  change  prediction  functions  due  to  air  drag 
with  a  standard,  baseline  ballistic  coefficient.  The  database  would  most  likely 
provide  coordinate  updates  in  terms  of  linear  functions  defined  by  intercepts 
and  slopes  for  various  classes  of  orbits.  Orbit  classes  should  be  determined  by 
a  combination  of  semi-major  axis,  eccentricity,  and  inclination.  Potentially,  the 
orbit  classes  could  also  be  defined  by  torus  coordinate  basis  frequencies. 

•  Study  the  relationship  between  ballistic  coefficient  vaiues  and  the  slope  of  the 
coordinate  time  derivative  prediction.  Based  on  the  author’s  brief  survey,  the 
relationship  appears  mostly  linear  with  small  quadratic  correction  terms. 

•  Develop  an  approach  for  “jumping”  to  an  adjacent  torus  when  momenta  differ¬ 
ences  become  too  large. 
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•  Incorporation  of  the  Moon’s  third-body  effects  as  additional  coordinate-momentum 
pairs  in  KAM  Torus  construction  should  be  investigated.  The  periodic  nature 
of  the  phase  accumulation  and  the  momenta  differences  seen  in  the  lunar  cases 
leads  the  author  to  believe  this  approach  is  promising.  Incorporating  the  Moon 
in  KAM  Torus  construction  would  also  remove  the  necessity  of  the  complex  op¬ 
erations  required  to  predict  the  Moon’s  state  over  time  for  input  into  third-body 
force  calculations. 
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Appendix  A.  Additional  Figures  and  Tables 


Table  11:  Comparison  of  Portion  of  Jacobian  Approximated  by  Two-Body  Problem 
and  Exact  Values  from  Torus  #2 


Two-Body  Jacobian 

Partial  Torus  Extraction  Jacobian 

3.672096e-l 

8.270228e-2 

3.661185e-l 

3.674205e-l 

8.271854e-2 

3.668758e-l 

-9.654751e-l 

-6.826067e-l 

-9.638614e-l 

-9.664381e-l 

-6.824581e-l 

-9.655660e-l 

2.132824e-l 

5.787414e-14 

2.134467e-l 

2.133621e-l 

2. 30448 le- 5 

2.133971e-l 

6.326216e-l 

8.924565e-l 

6.318474e-l 

6.319386e-l 

8.924208e-l 

6.315438e-l 

7.664626e-2 

3.394376e-l 

7.717644e-2 

7.678464e-2 

3.394487e-l 

7.700849e-2 

-7.390048e-l 

-3.893083e-14 

-7.384972e-l 

-7.394273e-l 

1.266063e-4 

-7.391464e-l 
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Figure  41:  Time  Derivative  of  Unperturbed  Torus  Coordinate  Utilizing  the  Two- 
Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  42:  Time  Derivative  of  Unperturbed  Torus  Coordinate  Q 3  Utilizing  the  Two- 
Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  43:  Time  Derivative  of  Unperturbed  Torus  Momentum  P2  Utilizing  the  Two- 
Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  44:  Time  Derivative  of  Unperturbed  Torus  Momentum  P3  Utilizing  the  Two- 
Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  45:  Time  Derivative  of  Unperturbed  Torus  #2  (ISS)  Coordinate  Q2  Utilizing 
the  Two-Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  46:  Time  Derivative  of  Unperturbed  Torus  #2  (ISS)  Coordinate  Qs  Utilizing 
the  Two-Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  47:  Time  Derivative  of  Unperturbed  Torus  #2  (ISS)  Momentum  Pi  Utilizing 
the  Two-Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  48:  Time  Derivative  of  Unperturbed  Torus  ^2  (ISS)  Momentum  P2  Utilizing 
the  Two-Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 


z> 


0.00004 

0.00003 

0.00002 

0.00001 

0.00000 

-0.00001 

-0.00002 

-0.00003 

-0.00004 


Time  [TU] 


Figure  49:  Time  Derivative  of  Unperturbed  Torus  #2  (ISS)  Momentum  P3  Utilizing 
the  Two-Body  Jacobian  and  the  Partial  Torus  Extraction  Jacobian 
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Figure  50:  Coordinate  Q2  Phase  Differences  for  Baseline,  Drag,  and  Lunar  Cases 


-Q 

0.00006 

tv 

0.00004 

a 

■a. 

0.00002 

'w' 

0.00000 

<U 

_c 

-0.00002 

IB 

-0.00004 

00 

-0.00006 

0.00006 

tT 

0.00004 

■ — 1 

0.00002 

Or 

e 

’O 

0.00000 

W) 

-0.00002 

Q 

-0.00004 

-0.00006 

0.00008 

tT 

0.00006 

0.00004 

Sr 

0.00002 

g 

0.00000 

-0.00002 

IV 

c 

-0.00004 

3 

—1 

-0.00006 

-0.00008 

Time  [days] 


Figure  51:  Coordinate  Phase  Differences  for  Baseline,  Drag,  and  Lunar  Cases 
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Figure  52:  Momentum  Differences,  8P2,  for  Baseline,  Drag,  and  Lunar  Cases 
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Figure  53:  Momentum  Differences,  5P3,  for  Baseline,  Drag,  and  Lunar  Cases 
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Figure  54:  Coordinate  Phase  Per  Time  and  Linear  Fits  for  280  days 


Figure  55:  Momenta  Differences  and  Linear  Fits  for  280  Days 
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Figure  56:  Coordinate  Qs  Phase  Differences  and  Phase  History  for  Lunar  Case 
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Figure  57:  Momentum  Differences,  8P\,  for  Baseline  and  Lunar  Cases 
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Figure  58:  Momentum  Differences,  SPz,  for  Baseline  and  Lunar  Cases 


-0.000003, 


0.00005 

0.00000 

-0.00005 

-0.00010 

-0.00015 

-0.00020 

-0.00025 

-0.00030 

-0.00035 

-0.000400L 


Time  [days] 


Figure  59:  Momentum  Differences,  SP>,  for  Baseline  and  Lunar  Cases 
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Figure  60:  Total  Position  Error  for  Reference  Torus  and  Perturbed  Motion  for  Base¬ 
line  and  Lunar  Cases 
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Figure  61:  Total  Velocity  Error  for  Reference  Torus  and  Perturbed  Motion  for  Baseline 
and  Lunar  Cases 
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